Kolmogorov algorithm for isochronous Hamiltonian systems

We present a Kolmogorov-like algorithm for the computation of a normal form in the neighborhood of an invariant torus in `isochronous' Hamiltonian systems, i.e., systems with Hamiltonians of the form $\mathcal{H}=\mathcal{H}_0+\varepsilon \mathcal{H}_1$ where $\mathcal{H}_0$ is the Hamiltonian of $N$ linear oscillators, and $\mathcal{H}_1$ is expandable as a polynomial series in the oscillators' canonical variables. This method can be regarded as a normal form analogue of a corresponding Lindstedt method for coupled oscillators. We comment on the possible use of the Lindstedt method itself under two distinct schemes, i.e., one producing series analogous to those of the Birkhoff normal form scheme, and another, analogous to the Kolomogorov normal form scheme in which we fix in advance the frequency of the torus.


Introduction and background
In various contexts in the literature, the use of the term 'Lindstedt series' for isochronous Hamiltonian systems often refers to one of two distinct methods, both applicable to the perturbative study of the dynamics around systems with elliptic equilibria. The difference between these two methods can be conveniently explained with the help of the following example: consider a 'Henon-Heiles' type Hamiltonian H = 1 2 (p 2 x + p 2 y ) + 1 2 ω 2 0,1 x 2 + 1 2 ω 2 0,2 y 2 + εP 3 (x, y) where, contrary to the actual Hénon-Heiles model ( [9]) (where ω 0,1 = ω 0,2 = 1), we first assume that the frequencies (ω 0,1 , ω 0,2 ) satisfy no resonance condition. P 3 can be any polynomial cubic in x, y.
The quantities ω i,1 (A 1 , A 2 ), ω i,2 (A 1 , A 2 ) are functions depending on two parameters A 1 , A 2 , called hereafter the 'amplitudes' of the oscillations in x and y respectively. They enter into the calculation through the choice made for the zero-th order terms x 0 (t), y 0 (t), since the iterative procedure starts by setting x 0 (t) = A 1 cos(ω 1 t + φ x0 ), y 0 (t) = A 2 cos(ω 2 t + φ y0 ) (4) where the initial phases φ x0 , φ y0 can be arbitrary. The above are common elements of the point of departure for both versions discussed below of the Lindstedt method. However, at this stage emerges an important bifurcation in the way we define the iterative scheme by which the functions x i (t), y i (t), ω i,1 , ω i,2 are to be computed. We discuss two distinct possibilities, referred to below as (i) a Lindstedt scheme 'analogous to the Birkhoff series', or (ii) a Lindstedt scheme 'analogous to the Kolmogorov series'.
As typical in perturbation theory, the formal difference between the above two schemes actually reflects a real (physical) difference in the way we interpret the meaning of the series (3). In summary, the difference can be posed as follows (see section 2 for details): (i) in the scheme called below 'analogous to Birkhoff', we seek to construct a quasi-periodic solution valid for any value of the amplitudes A 1 , A 2 within a suitably defined open domain around the origin. Thus, the series (3) in this scheme are meant to answer the question of what are the values of the frequencies ω 1 , ω 2 under which the motion takes place for any given and pre-selected sets of values of the amplitudes A 1 , A 1 in the above domain. The reader is referred to [15] where a clear exposition of the method is given in the framework of special solutions of the three-body problem computed via Lindstedt series.
(ii) in the scheme called below 'analogous to Kolmogorov', instead, we fix in advance the values of the frequencies ω 1 , ω 2 (see [13] for a clear exposition of the method in the context of the forced anharmonic oscillator); this is called by some authors a 'torus fixing method'. A relevant remark in the context of this last method is that the series (3) are actually purported to answer the question reverse to the one posed in (i) above. That is, the question now is: with given and pre-selected values of the frequencies ω 1 , ω 2 , invert the series (3) and compute which are the corresponding amplitudes A 1 , A 2 for which we obtain quasi-periodic trajectories with the frequencies ω 1 , ω 2 . Thus, in method (i) the series are parameterized by the amplitudes A 1 , A 2 , which can be selected at the beginning of the construction, while in method (ii) the solutions are parameterized by the frequencies ω 1 , ω 2 , which are the parameters to select at the beginning of the construction. Also, in the latter case the series inverse to (3) turn out to have the form (in the cubic case) for some constant coefficients C i,1 , C i,2 computable from the series (3). Thus, with all frequencies of the problem fixed in advance, establishing the convergence of the inverse series (5) suffices to answer the question posed at (ii). The question of the convergence of the series is, of course, crucial, and related to the kind, and pattern of accumulation in the series terms, of small divisors appearing at successive perturbative steps. As regards the kind of divisors, we can readily see that: -in scheme (i) we obtain divisors of the form k 1 ω 0,1 +k 2 ω 0,2 , with (k 1 , k 2 ) ∈ Z 2 , |k 1 |+|k 2 | = 0. This follows from the kind of linear (non-homogeneous) equation to solve iteratively. Deferring details to the example treated in Section 2, we briefly recall that in scheme (i) we introduce the parametrization (modulo two unimportant phases) ϕ 1 = ω 1 t, ϕ 2 = ω 2 t, and after introducing the series expressions (2) and (3) to the equations of motion and separate terms of like orders we arrive at equations (to be solved iteratively) of the form: where the functions Θ 1,i (ϕ 1 , ϕ 2 ), Θ 2,i (ϕ 1 , ϕ 2 ) contain trigonometric terms in the angles ϕ 1 , ϕ 2 (see [15], section 4).
-In scheme (ii), instead, we obtain divisors of the form k 1 ω 1 + k 2 ω 2 , i.e., depending on the (fixed) pre-selected new frequencies ω 1 , ω 2 . This follows from the fact that the linear non-homogeneous equations to solve are now of the form (see [13]): containing trigonometric terms in the angles ϕ 1 = ω 1 t, ϕ 2 = ω 2 t. Note that since the divisors depend on the new frequencies ω 1 , ω 2 , choosing non-resonant values for the latter permits the formal construction to proceed; this, even when the unperturbed frequencies ω 0,1 , ω 0,2 are, instead, resonant. As regards convergence, in the case (i) Poincaré ([19], Ch.IX) already emphasizes that the Lindstedt series with divisors depending on the original harmonic frequencies ω 0,1 , ω 0,2 are divergent, exhibiting the well known asymptotic character associated with the series computed via a Birkhoff normal form (see [5] for a review). Indeed, as shown by example in section 2 below, it possible to construct Birkhoff series yielding the same individual solutions as those of the Lindstedt series of scheme (i). We note here that the series originally introduced by Lindstedt ( [16,17,18]), albeit somewhat different in structure, exhibit the same divisors as those of the scheme (i) above, thus, according to Poincaré, they are only asymptotic. On the other hand, Eliasson ([6]) and Gallavotti ([11,12]) established the existence of convergent Lindstedt series by the 'torus fixing method' on the basis of the cancellations between terms with small divisors (see [14] for an instructive example). A proof of the convergence of scheme (ii) is actually possible by diagrammatic methods via the following theorem [4]: 4]). Consider the N coupled oscillator equations where ε is a real parameter, f (x, ε) is real analytic at x = 0 , ε = 0 and at least quadratic in x, and the frequency vector ω is diophantine. Let be a solution of (8) for any choice of the complex constants c j and for ε = 0. Let Γ(c) = max(|c 1 |, . . . , |c N |, 1). Then, there exists a positive constant η 0 and a function η(ε, c) holomorphic in the domain |ε|Γ 3 (c) ≤ η 0 , real for real ε, such that the system admits a solution of the form holomorphic in the domain |ε|Γ 3 (c)e 3|ω|| t| ≤ η 0 and real for real ε, t. The constants A ν are O(ε), except for the constants A 1,0,...,0 , A 0,1,...,0 , A 0,0,...,1 , which are equal to c 1 , c 2 , . . . , c N respectively.
A similar proof in action-angle variables in the case N = 2 is discussed in [1]. As a final introductory remark, the series construction in the isochronous case finds a plethora of applications in various fields of physics. We mention in particular, the use of the Lindstedt method for the computation of solutions lying on low-dimensional tori ('q-tori') in the celebrated Fermi-Pasta-Ulam (FPU) problem ( [2,3]). The FPU model takes, in normal mode space, the form of N harmonic oscillators coupled with nonlinear terms:Q where the frequencies Ω k , k = 1, . . . , N are given in terms of the FPU normal mode spectrum Ω k = 2 sin(kπ/(2(N + 1)), the function F can be cubic or quartic in the variables Q k , and the perturbation ε satisfies some scaling law with N . Flach and co-workers ( [7,8]) emphasized the special role for dynamics played by solutions called 'q-breathers'. These are periodic orbits of the form Q q (t) = A q cos(ω q t + φ q ), and Q k (t) = 0 for k = q. For the frequencies ω q we obtain series expressions of the form Then, for ε sufficiently small, the Lindstedt method (ii) above allows to represent the q-breathers via the Fourier expansion wheref k,m = O(ε p(k,q,m) ), with integer exponent p(k, q, m) ≥ 1. The relevant point for the FPU problem is that the rules of propagation of the amplitude A q in the series terms for all modes allows to find an analytic formula explaining the phenomenon of 'energy localization' observed for particular initial excitations in the FPU model. In [2] and [3], on the other hand, it was shown that the q-breathers constitute only the first member in the hierarchy of special FPU solutions that exhibit energy localization. More general members are the 'q-tori', i.e., special solutions with M < N incommensurable frequencies satisfying where R q = (q 1 , . . . q M ) ∈ {1, 2, ..., N } M . The corresponding Fourier representation of these special solutions can again been computed using Lindstedt series, and it obtains the form . . , φ q M ), and f k,m = O(ε) for all k = 1, . . . N . Furthermore, the propagation of the amplitudes A k in the series terms allows to interpret a variety of complex localization profiles encountered for particular initial mode excitations in the FPU problem (see the corresponding theorems in [3]). We mentioned already that for the Lindstedt series analogous to the Birkhoff ones there exists a Birkhoff normal form yielding the same solutions as those recovered by the Lindstedt method via an indirect approach, i.e., one based on a sequence of normalizing transformations involving canonical changes of variables. It is natural to ask whether this correspondence between a direct (Lindstedt) and indirect (normal form) method extends in the case of the torus-fixing method as well. Due to the lack of a twist condition, the torus-fixing process in the isochronous case has to be dealt with using a technique based on 'counterterms' (see [11]), or a KAM algorithm 'with knobs' (see [20] in the present volume). We have developed a Kolmogorov algorithm using counterterms, which is able to recover the solutions of the direct Lindstedt method in both cases of full or low-dimensional tori. In the present paper we will discuss in some detail the proposed Kolmogorov algorithm for the fulldimensional case only, deferring the low-dimensional case to a different publication. Section 2 gives an elementary example allowing to fix all principle ideas in this comparison. Section 3 gives the general algorithm for the Kolmogorov method in the isochronous case. This is similar as the algorithm 'with knobs' discussed by Sansottera and Danesi ([20]) in the present volume. We point out, however, some differences in our viewpoints regarding what the algorithm actually achieves to compute and how it should be implemented in the context of isochronous models with given initial unperturbed frequencies ω 0 . Some results and detailed proofs are deferred to the Appendix.

An elementary example
In order to illustrate the methods discussed above, we will consider an elementary example stemming from the following one-degree of freedom Hamiltonian with a even power dependence on the canonical variables Using the harmonic oscillator action-angle variables (J, q) with x = √ 2J sin(q) , p = √ 2J cos(q) , we obtain Let J 0 be the label of a given torus (periodic orbit) of the harmonic oscillator model H 0 . Consider a real neighborhood D ε = {J = J 0 + p with |p| < D ε }, where D ε = O(ε). We will illustrate four different perturbative methods to treat the dynamics in the phase-space neighborhood T × D ε : these are i) a Birkhoff normal form construction, with ii) its analog in terms of Lindstedt series, iii) a Lindstedt series exhibiting the torus-fixing property of the Kolmogorov method, and, finally iv) the normal form analogue of iii). In the next section, we will give the general formulas for the construction of the Kolmogorov normal form in the case of n-degree of freedom Hamiltonians with isochronous integrable part H 0 .

Birkhoff normal form
Setting J = J 0 + p the Hamiltonian takes the form (apart from a constant) A 'Birkhoff normal form' for the Hamiltonian (19) can be computed by introducing a canonical transformation eliminating the angle q . Writing the Hamiltonian (19) as where Z 0 = ω 0 p and f 1 = k,l c k,l p k e ilq we will define a Lie generating function χ (1) (q, p) bringing this Hamiltonian to normal form up to terms O(ε). In the standard procedure, it is sufficient to set χ (1) = X (1) , where X (1) solves the homological equation L X (1) Z 0 + εf 1 = εζ 1 , with ζ 1 =< f 1 > q and L X denoting the Poisson bracket operator L X = {·, X}. However, comparing the result with the one obtained by the corresponding Lindstedt method (see next subsection) requires a small modification in the definition of χ (1) . Consider the canonical transformation (q, p) → (q (r) , p (r) ) obtained after r normalization steps: we require that the canonical transformation be such that the initial condition q (r) = p (r) = 0 in the new variables be mapped to the initial condition q(q (r) = 0, p (r) = 0), p(q (r) = 0, p (r) = 0) = (0, 0) + O(ε r+1 ) in the original variables. It is easy to see that such a requirement of control on the initial condition can be fulfilled by setting where K (n) , S (n) are constants possible to compute at every step by requiring that Since {ω 0 p, K (n) q} = −ω 0 K (n) , this procedure will only alter the normal form at the n-th step by a constant, adding, however, some trigonometric terms to the remainder at every step. As an example, we can readily verify the following formulas for the first step (and analogously for subsequent steps) Omitting details, the formulas obtained after two normalization steps as above are the following: the Hamiltonian takes the form (2) p=p (2) . For simplicity in the notation, from now on we omit superscripts from the variables q and p unless explicitly required, adopting, instead, the convention that the symbols (q,p) in any function of the form F (r) (q,p) imply the new canonical variables computed after r normalization steps. Then, up to order 2 in ε we obtain: where the remainder R (2) is O(ε 3 ). The Hamiltonian Z (2) can now be used to analytically compute ε 2 −precise solutions to the equations of motion in the variables (q,p). The equations of motion are fixing the initial conditionp(0) =q(0) = 0 yields the solutioñ This can be back-transformed to the solution in the original variables. We have Substituting the solutions (q(t),p(t)) in the previous expression, we find Observe that q(0) = 0 and J(0) = J 0 , as was required. Two remarks are in order: (i) The divisors appearing in all series expressions obtained above depend on the (unique, in the case of 1DOF systems) unperturbed frequency of the model, i.e., the frequency ω 0 of the linear oscillator. In the case of systems with N > 1 degrees of freedom, the above series will produce divisors of the form m · ω 0 , where m ∈ Z N , |m| = 0, and ω 0 is the N −vector of the unperturbed frequencies of the problem. This implies that the method may formally proceed only when the vector ω 0 is non-resonant. As regards the series convergence, this is guaranteed in an open domain in the 1DOF case. However, when N > 1 the series are in general only asymptotic (see [5] for reviews).
(ii) The value of the initial datum for the action J 0 can be chosen at the end of the process, i.e., the numerical value of J 0 need not be fixed in advance. In fact, J 0 (or, more generally, the vector J 0 ) can be regarded as a parameter free to transfer in all iterative computations. It is, then, in terms of this parameter that we obtain an expression for the frequency(ies) ω on the torus with the initial conditions J(0) = J 0 , q(0) = 0 as a function of J 0 . In our example this function is given up to O(ε 2 ) in Eq (20).
(iii) Any choice is possible to impose on the initial phase q(0), altering slightly the definition of the constants S (n) .

Lindstedt solution analogous to Birkhoff
Analytical solutions in the form of formal series obtained like (22) above will be hereafter called 'Birkhoff solutions'. This terminology emphasizes the fact that the solutions are obtained using a Birkhoff normal form procedure. It is notworthy that Poincaré's reference to the 'Lindstedt series' (Methodes nouvelles, Ch. IX) actually consists of the construction of solutions to coupled oscillator problems obtained 'indirectly', i.e., via transformations of the variables and the engineering of the corresponding Hamiltonian function. Hence, the more general term 'Poincaré-Lindstedt' series often encountered in literature. We now give, instead, the 'direct' (i.e., without transformations) series construction of the Birkhoff solution (22) by implementing, instead, the original method of Lindstedt in the framework of the canonical action-angle variables of the harmonic oscillator model.
The Hamilton's equations of motion for the Hamiltonian (18) are We perform the following steps: Step 1: time re-parametrization. Set ϕ = ω t where ω is the (still unknown) frequency on the 1-torus (periodic orbit) corresponding to the solution with initial conditions J(0) = J 0 , q(0) = 0.
Step 2: frequency expansion. This is the key element of the Lindstedt method. We expand ω as Note that the corrections a i i = 1, 2, . . . are functions of the parameter J 0 , i.e., of the 'amplitude' (square) of the oscillations.
Step 3: expansion of the solution. We write Note here a key element of the method, which is the fact that all functions q i , J i , i = 0, 1, . . . are considered functions of the phase ϕ = ω t rather than of the time t itself. This is the key point in the differentiation of the method presented here with respect to the 'torus-fixing' method presented in Subsection 2.3. The relevant fact is that the use of the angle ϕ (or, more generally, of a set of angles ϕ ∈ T N , with ϕ = ωt in the N-DOF case), instead of time, allows to split the equations of motion in a sequence of linear non-homogeneous equations, whose iterative solution introduces divisors depending on the unperturbed frequencies ω 0 . This is made clear in the next step.
Step 4: splitting of the equations of motion in powers of ε and iterative solution. In our example, from the definition of ϕ we have: Substituting (24) and (25) in the equations of motion (23a) and (23b) leads to the following expressions (up to order 2 in ε ) Now, in order to iteratively determine the functions q i (ϕ), J i (ϕ), we compare the two sides of the equations of motion at equal orders. At order 0 we have Fixing the initial data as q 0 (0) = 0 and J(0) = J 0 , we then arrive at q 0 (ϕ) = ϕ and J 0 (ϕ) = J 0 . At order 1, we now have As well known, a 1 is determined by the requirement that no secular terms be present in the series. This yields a 1 = −3 J 0 /4 . Then, the Cauchy problem with the previous differential equations and the initial constants q 1 (0) = J 1 (0) = 0 yields the solution The process can be repeated at subsequent orders. We leave to the reader to verify the result at second order, leading eventually to Comparing the result with Eqs (20) and (22), we can see that the solutions q(t) , J(t) of the two methods are equal. This is easy to justify by checking the structure of the l.h.s of Eq (26). After the expansion of the frequency, the differential equations to solve at subsequent steps all involve the operators ω 0 d/dϕ. Hence, all divisors appearing in the series terms are in terms of the unperturbed frequency ω 0 .

Lindstedt series analogous to Kolmogorov
It was already pointed out that the Lindstedt series examined so far, as well as their 'indirect' (normal form) counterpart, produce, in general, series which are divergent and only have an asymptotic character. 2 On the other hand, Eliasson [6] and Gallavotti [12] established the existence of convergent Lindstedt series in nonlinear Hamiltonian systems satisfying the necessary conditions for the holding of the Kolmogorov-Arnold-Moser theorem. Gallavotti ([11]) presented a diagrammatic proof of the convergence of the Kolmogorov normal form also in the 'twistless' case, i.e., when the size of the Hessian matrix of the unperturbed Hamiltonian H 0 (J) is not limited from below (it is equal to zero in the 'isochronous' case). A constructive Kolmogorov algorithm able to deal also with the twistless case is presented, along with the demonstration of its convergence, in the present volume by Sansottera and Danesi [20]). The convergence of the Lindstedt series in this case is addressed, instead, in [4] (see remarks in the introduction). As discussed in the introduction, when series constructions analogous to Kolmogorov are sought for in the isochronous case, an important point to address is the need for performing, at the final stage of the construction, a process involving series reversal. This reversal is necessary in order to explicitely compute the solutions q(t), J(t) whose initial conditions q 0 , J 0 correspond to motion on a torus with given frequency vector ω. The key remark is that the value of J 0 , which parametrises the solutions, cannot be fixed in advance due to the lack of a twist condition allowing to compute the mapping ω(J 0 ). We now discuss the application of the direct (Lindstedt) method analogous to Kolmogorov in the same example as in the previous two sections, aiming to illustrate the above points.
Fixing the frequency of the torus in the isochronous case can be implemented as described in [13]: assume that we target a particular solution of the equations of motion (23a) and (23b) represented as a trigonometric series and evolving according to a given pre-selected frequency ω. Inverting the expansion (24) we obtain Also, as before we expand the solution as Note, however, that this time we perform no time-reparametrization, i.e., the solutions remain expressed as functions of the time t. Thus, replacing the above expressions in (23a) and (23b) the equations of motion lead now to the expressions (up to order 2 in ε ) cos(4 q 0 (t)) .
Collecting now the terms of equal order we can compute iteratively all the functions q i (t), J i (t), for i = 0, 1, 2, setting, at order zero: q 0 (t) = ω t , J 0 (t) = J 0 which corresponds to the choice of the initial condition q 0 (0) = 0 and J 0 (0) = J 0 . Then, at first order we haveq implying a 1 = −3 J 0 /4 and Note that the integration constants in (32) were set as q 1 (0) = J 1 (0) = 0 , consistent with our choice of initial condition. Repeating the procedure at second order yields It is instructive to compare the solutions (33) above with those (Eq (28)) found in the case of the original Lindstedt method. We have the following remarks: (i) From second order and beyond the coefficients in front of the harmonics cos(mωt), sin(mωt) for the same m are not all equal in the two solutions.
(ii) Most importantly, the structure of the divisors in the two solutions is different. In the torus fixing case (Eq (33)), all divisors involve the corrected frequency, ω, instead of the original frequency ω 0 of the linear oscillator. This is not a problem for the method to proceed, since this frequency is known in advance. In N −DOF systems, in general, with the 'torus fixing method' we obtain divisors of the form m · ω. Thus, the method can proceed even when the original frequencies of the N linear oscillators ω 0 are resonant, as long as we impose a non-resonant detuning in the adopted values for the final frequencies ω on the torus.
(iii) On the other hand, the value of the initial condition J 0 leading to motion on the torus with frequency ω remains uknown up to the end of the construction. From the point of view of the symbolic implementation of the method in the computer, J 0 is a symbol whole powers have to be carried on along with the remaining powers of trigonometric monomials in all series terms and at all iterative steps.
(iv) At the end of the process, however, J 0 can be estimated by reversing the series (29), given that the functions a n (J 0 ) are monomials of degree n in J 0 . In an analogous way, in the N −DOF case we will end up with N series equations of the form ε n a n,j (J 0,1 , . . . , J 0,N ) where the functions a n,j (J 0,1 , . . . , J 0,N ) have all been specified iteratively up to a maximum truncation order in n, and they are polynomial in J 0 ≡ (J 0,1 , . . . , J 0,N ). Thus, the series (34) can be formally inverted, yielding where the functions P n,j (ω − ω 0 ) are polynomial of degree n in the differences (ω − ω 0 ). Note that for the inverse series to converge we must require the difference |(ω − ω 0 )| to be smaller than ε, a fact which limits how far we can detune ω from ω 0 in order to be able to specify the corresponding initial condition J 0 . An obvious choice is |(ω − ω 0 )| = O(ε 2 ). At any rate, we emphasize that the convergence of the inverse series (34) is an open problem, crucial to the applications. 3

Kolmogorov normal form
We finally arrive at the here proposed Kolmogorov normal form algorithm yielding solutions equivalent to those discussed in the previous subsection. This is implemented by the following steps: Step 1: substitution of the frequency series into the Hamiltonian. In our example, we substitute ω 0 in the Hamiltonian (19) with the series (29). This yields (up to second order, apart from constants) Step 2: normalization. To set the Hamiltonian into Kolmogorov normal form up to second order, we fix the value of the constants a 1 , a 2 and we perform a sequence of Lie transformations aiming to give the Hamiltonian the form (in the transformed variables) To this end: -First order: we fix a 1 so that the linear term εa 1 p acts as counterterm for the term ε(3J 0 /4)p. This provides, in the twistless case, a process by which the frequency ω can be kept fixed (in the usual twist case, instead, this would have been accomplished by exploiting the Hessian matrix of H 0 ). Formally, we require that leading to a 1 = − 3 J 0 4 . Now, we insert this expression for a 1 in the Hamiltonian (36), leading to We can now eliminate the O(ε) trigonometric terms in the Hamiltonian with the usual procedure. Namely, we define a generating function X (1) (q) used to eliminate terms not depending on the action p . These are Note that, similarly as in the Birkhoff case (Subsection 2.1), here too we have to fix the initial conditions so that the relation q(0) = p(0) = 0 is preserved between variables before and after the canonical transformation. This is achieved by setting the generating function as χ is a constant. The general rules for the determination of the constants K (j) will be discussed in Section 3.
Using the generating function χ 1 (q) we obtain the intermediate Hamiltonian We will now eliminate the trigonometric terms linear in the momentum in H (1) . These are 1,1 = 0 . As before, the value of the constant S (1) is fixed by the requirement that the initial phase of the solution of q be preserved to zero by the transformation. We find S (1) = 0 (all constants S (i) = 0 when the Hamiltonian has an even symmetry). Theñ

The new Hamiltonian is
and it is in Kolmogorov normal form up to order ε .
Repeating the same procedure at second order we arrive at the formulas leading to the Hamiltonian H (2) (q,p) = Z (2) (q,p) + R (2) (q,p) with the Kolmogorov normal form part 64 ω cos(8q) ((q,p) denote again the new variables after new normalization steps).
Step 3: calculation of the solution on the torus. Using the compact notation where R j (q,p) = O(||p || 2 ) j = 1, 2 , the equations of motions under the Hamiltonian Z (2) are .
The torusp(t) = 0 ,q(t) = ω t (where we choseq(0) = 0 ) is a solution of this system. This can be back-transformed in the original variables using Substituting the solution (q(t) = ω t,p(t) = 0) in the previous expressions, we readily deduce the solution in the original variables: 32 ω 2 sin(2 ω t) + (42) Also, using the computed expressions of a 1 and a 2 (Eqs (37) and (40)), we obtain the relation between the torus frequency and the original frequency These expressions are identical to the ones found by the Lindstedt method of Subsection 2.3 (see (33)).

Comparisons and numerical tests
In order to better visualize the differences between the methods (i) and (ii) (i.e., respectively, 'analogous to Birkhoff' and 'analogous to Kolmogorov') we report, in the Tables 1 and 2 below, the series terms corresponding to the solutions for q(t) and J(t) as obtained by the two methods up to order ε 2 . Table 1: Comparison between the series terms for the solution q(t) in the Lindstedt series obtained by the methods (i) and (ii) up to order ε 2 .
Method:  Method:  Table 1 above, we observe that the solutions q(t) and J(t) obtained by methods (i) and (ii) have the same form up to order O(ε) , except for the fact that in method (i) all divisors depend on the frequency ω 0 rather than ω . On the other hand, the coefficients found by the two methods start differing from the order ε 2 and beyond. Moreover, the series for the frequency ω obtained by the two methods are: which differ again as regards their divisors. The effect of these differences on the precision of the two methods can be seen even in the simplest case of a system with one degree of freedom in which both series are convergent (the series 'analogous to Birkhoff' are, instead, divergent when N ≥ 2 ). To this end, we report below a comparison between the numerical solutions and the analytical ones, obtained by the methods (i) and (ii) in the example of the Hamiltonian (17). To produce the Lindstedt series (ii) in this example we work as follows: fixing the initial and final frequency ω 0 and ω , we reverse the series (24) according to where dω = ω − ω 0 and f −1 denotes the series inverse to f . Then, having specified J 0 through the inverse series (43), we compute all numerical coefficients in the Lindstedt series (i) and (ii) up to order ε 4 . We analyze three different cases: 1) ε = 1 , ω 0 = 1 and ω = 1.002 ; 2) ε = 1 , ω 0 = 1 and ω = 1.02 ; 3) ε = 1 , ω 0 = 1 and ω = 1.2 .
Let us note that all results are rescalable to different choices of ε . Reversing the series (24) as prescribed by (43), we obtain, respectively, the following values for the amplitude (initial action datum): J 0 = 0.0026769 , 0.0277048 , 0.383509 . Substituting these values of ε , ω , ω 0 and J 0 in the solutions q(t) and J(t) we obtain the solutions, as functions of t , produced by methods (i) and (ii). At the same time, it is possible to integrate the equations of motion (23a) and (23b) to produce a numerical solution, starting from the initial conditions q(0) = 0 , J(0) = J 0 . Figure 1 shows the difference (in Log 10 scale) between the semi-analytical solutions obtained by methods (i) and (ii) as above and the numerical solution for the angle q(t) . Since the only systematically growing error is due to differences in the frequency estimates, all errors between the numerical and analytical solutions in the angle q(t) grow linearly in time. We observe, however, that the Lindstedt method 'analogous to Kolmogorov' always produces more precise results than the one 'analogous to Birkhoff' with a difference in precision of about one order of magnitude when dω is of order 10 −3 and raising up to two orders of magnitude when dω becomes of the order of unity.

Hamiltonian preparation in the case of odd nonlinear couplings
A particularity of the example treated above is the fact that the original Hamiltonian is analytic in the whole domain J ∈ R. This changes, however, in more general models in which the Hamiltonian is of the form H(q, J ) = ω 0 · J + εh(q, J ; ε) where the development in powers of the variables J contains semi-integer powers, as is, for example, the case of polynomial nonlinear couplings containing odd terms in one or more of the oscillator variables Figure 1: Difference between the numerical solution q(t) and the one produced by the method 'analogous to Birkhoff' (in blue, (i)) or 'analogous to Kolmogorov' (in red, (ii)). From left to right for the cases 1) ε = 1 , ω 0 = 1 and ω = 1.002 , 2) ε = 1 , ω 0 = 1 and ω = 1.02 , 3) ε = 1 , ω 0 = 1 and ω = 1.2 .
x j , y j . Let us note that physical examples of such Hamiltonian systems, with oscillators non-linearly coupled through odd polynomial terms, are ubiquitous, and include the Fermi-Pasta-Ulam α−model, the secular Hamiltonian in resonant cases of perturbed Keplerian N −body problems, magnetic bottle Hamiltonian models, etc. In the next section, we will present a formal algorithm applicable to a generic form for the function H 1 . The sequence of normalizations in this algorithm is arranged so that the result agrees with the corresponding one obtained with the Lindstedt series. A particular example showing this agreement is presented in the Appendix.

KAM algorithm for isochronous systems
We now give the Kolmogorov algorithm for generic isochronous systems with Hamiltonian where q ∈ T n , J ∈ B ⊆ R n ,h i, s = O(||J || s 2 ) . We assume that the Hamiltonian (44) has a half-integer power dependence on J , i.e. admits the expansion, for n → ∞, of the truncated series whereh (k) are the k−th derivatives ofh with respect to p = J − J 0 . Apart from constants, the Hamiltonian can be written in the compact notation The algorithm allows to compute quasi-periodic orbits with a frequency ω fixed in advance, given by where the parameters J 0 , whose values are to be specified in the end of the process, give the initial conditions for J of a trajectory on the torus with frequencies ω . To this end, along with the normalizing canonical transformation, the algorith computes 'on the go' the functions a i (J 0 ) ('counter-terms'). Substituting the series (47) in the Hamiltonian (46) we arrive at: In the following we use the notation where h (k) i,0 is independent from p , h  such that, after r normalization steps, the Hamiltonian (48) is given by the formal series: Given the Fourier series h Proof of the proposition. The generic r-th iterative step of the algorithm is defined as follows: after r − 1 steps, the hamiltonian (48) has the form: where C 1 . . . C r−1 are constants and h (r−1) i ∈ O(|| p || 2 ) ∀ i = 1, . . . , r − 1 . Taking into account the notation (50), we re-define the quantity h We then have: First part of the proof: the r-th normalization step consists of two substeps, each involving a canonical transformation.
• First half step: we set where the generating function χ (r) 1 is defined as χ (r) 1 (q) = X (r) (q) + K (r) · q , with K (r) an appropriate constant vector defined below, and X (r) is defined through the homological equation: where . denotes the mean over q . The function X (r) (q) eliminates all terms of order O(ε r ) depending only on the angles q in the Hamiltonian H (r−1) . Writing h (r−1) r,0 (q) in the Fourier form h (r−1) r,0 (q) = k∈Z n c (r−1) k e ik·q from (54) we find: We then impose the condition that the terms of order O(ε r ) linear in p have zero This specifies a r as a function of J 0 via Eq (52). Finally, we specify the constant vector K (r) . To this end, we impose the condition that the solution in p has the form k A k (cos(k · q) − 1) , so that, at time t = 0, we have p(0) = 0. Writing (55) as the poisson bracket of χ Hence i.e., Finally, we compute H (r) in (53) as: where h (r) • Second half step (of r-th): We compute with a generating function χ jk p j e ik·q the solution of (70) is We finally compute the constant vector S (r) by the condition q(0) = 0 at the time t = 0. By the poisson bracket we obtain The Hamiltonian H (r) (Eq. (69)) is: where h (r) Equation (79), using homological equation (70), can be written equivalently as: where, from (79) we have where we used Eqs (62), (63) and (64), respectively. Using (68) we have Equation (81) can now be written in the form ∀ k ≥ r + 1 , i ≥ 2 . The indices i, j, k in the first term of the above equation satisfy the relation where we have used the inequalities j ≤ (k − 1)/r − 1 and k − 1 ≤ k − 1 . Then Eq (88) ensures that the first term of (87) satisfies the definition (67). For the second term in (87), we have the following useful (also in the sequel) remark: we can write k = mr + f , where 0 ≤ f ≤ r − 1 and m ∈ N . Thus Then We can conclude that 1 ≤ k − (k − 1)/r r ≤ r , implying that the definition (68) holds for the second term in (87). We can then write (87) as 4 h (r) ∀ k ≥ 2r + 1 . We study separately the relations satisfied by the indices (j, k) of each of the three terms in the previous equality. Following (88), we have that . Then, the definition (64) holds for the first term. For the second and third terms, we have different formulas according to whether or not k is a multiple of r.
(i) First case: k is not a multiple of r , i.e., From (91) we have that 1 ≤ k − k−1 r r ≤ r − 1 and, consequently, r + 1 ≤ k − k−1 r r + r ≤ 2r − 1 . Then, the definitions (62) and (63) hold, respectively, for the third and second term of (93). Thus, we can write equation (93)  (ii) Second case: k is a multiple of r , i.e., k = mr with m ∈ N .
From (90) we now have that k − k−1 r r = r and, consequently, k − k−1 r r + r = 2r . Then, the definitions (62) and (64) hold, respectively, for the third and the second term of (93). Thus, we can write Eq (93) as ∀ k ≥ 2r + 1 , k = mr , m ∈ N . By the same argument, we can write Eq (78) as: ∀ k ≥ r + 1 , k = m r (m ∈ N) . In view of the inequalities (88) and (91) we then readily find that the definitions (67) and (65)  (ω · p) , ∀ k ≥ r + 1 , k = m r (m ∈ N) (where we used the definitions (67) and (66) for the first and second part of the sum). As before, we can then write the previous equation as Use of Eq (66) and the homological equation (70) then concludes the proof.

Concluding remarks
The focus of the present paper is twofold.
i) We emphasize the differences between two distinct methods by which Lindstedt series can be computed in nonlinearly coupled oscillator Hamiltonian models. In particular, we discussed a Lindstedt method called 'analogous to the Birkhoff series', and another called 'analogous to the Kolmogorov series'. In the first case, series expansions are defined in open domains in the action space around the origin, and the frequencies are represented as (series) polynomial functions of the action variables. In the second method, instead, the frequencies along any required torus solution must be fixed in advance ('torus fixing'), and the series allow (after reversion) to determine a posteriori the amplitudes (or values of the actions) for which the motion takes place with a given set of frequencies. We explain how this difference in the physical interpretation of what the series are meant to compute leads also to formal differences in the way the series terms are computed, as well as to real differences in the convergence and precision properties of the series.
ii) We propose and give the formal structure of an algorithm viewed as an analogue of the Kolmogorov scheme within the framework of the 'torus fixing' method, but properly dealing with the lack of the twist property in the Hamiltonian appearing in the kernel of the associated perturbative scheme. In particular, we demonstrate how the lack of the twist condition can be compensated by a suitable introduction of counter-terms (in the spirit of the procedures proposed already in [1,11,12,20]). In fact, we demonstrate how to compute such counter-terms in a way leading to precisely the same frequency corrections as those obtained by the Lindstedt method called 'analogous to Kolmogorov'.
Our formal algorithm given in Section 3 above is accompanied by particular examples offering some intuition into the subtleties of each examined method. On the other hand, as stressed in the introduction, the convergence of the Kolmogorov series computed by the algorithm presented in Section 3 can only be inferred by an indirect argument, namely, their equivalence with the Lindstedt series of Subsection 2.3. Thus, we point out the interesting open question of a direct proof of the convergence of our hereby presented KAM algorithm for isochronous systems (see also [20]).