Kolmogorov algorithm for isochronous Hamiltonian systems

Foreword

The motivation for this paper stemmed from an exchange of e-mails and discussions that the authors of refs. ^[2] and ^[3] (C. E. and Eleni Christodoulidi) had with Antonio Giorgilli, Sergej Flach and Giancarlo Benettin several years ago. The discussion was on the meaning of Lindstedt series as a method to construct invariant tori in perturbed oscillator systems like the celebrated problem of Fermi-Pasta-Ulam (see Galgani ^[10] in the present volume). We, as many others, are fond of recognizing the deep influence of Antonio's thought on our way of viewing the role of perturbative methods in the understanding of Hamiltonian dynamical systems. Since this is a paper to appear in the special issue in his honor, our presentation below is, in turn, strongly shaped by what we hope could be a framework of discussion interesting to our inspiring and beloved teacher A. Giorgilli.

1.   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

where, contrary to the actual Hénon-Heiles model (^[9]) (where $\omega_{0, 1} = \omega_{0, 2} = 1$), we first assume that the frequencies $(\omega_{0, 1}, \omega_{0, 2})$ satisfy no resonance condition. $P_3$ can be any polynomial cubic in $x, y$.

We are interested in constructing perturbative series solutions in the model (1.1) under the form:

where we adopt some periodic form for the functions $x_0(t), y_0(t)$ and compute iteratively all subsequent functions $x_i(t), y_i(t)$, $i = 1, 2, \ldots$. We further want to secure that the iterative scheme used to compute the functions $x_i(t), y_i(t)$ preserves the quasi-periodic character of the solutions, i.e., produces no secular terms (of the form $t\sin(\omega t)$, etc) for some frequencies $\omega$ obtained as discussed below. An elementary remark in this context is that the nonlinear coupling of the oscillators implies that quasi-periodic orbits in the above model are expected to evolve, in general, with frequencies $\omega_1, \omega_2$ different from those of the unperturbed oscillators, $\omega_{0, 1}, \omega_{0, 2}$. As it is well known (see, for example, ^[13]), recognition of this fact implies to introduce formal series also for the frequencies:

The quantities $\omega_{i, 1}(A_1, A_2), \omega_{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

where the initial phases $\phi_{x0}, \phi_{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), \omega_{i, 1}, \omega_{i, 2}$ are to be computed. We discuss two distinct possibilities, referred to below as (ⅰ) a Lindstedt scheme 'analogous to the Birkhoff series', or (ⅱ) 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 (1.3). In summary, the difference can be posed as follows (see section 2 for details):

(ⅰ) 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 (1.3) in this scheme are meant to answer the question of what are the values of the frequencies $\omega_1, \omega_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.

(ⅱ) in the scheme called below 'analogous to Kolmogorov', instead, we fix in advance the values of the frequencies $\omega_1, \omega_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 (1.3) are actually purported to answer the question reverse to the one posed in (ⅰ) above. That is, the question now is: with given and pre-selected values of the frequencies $\omega_1, \omega_2$, invert the series (1.3) and compute which are the corresponding amplitudes $A_1, A_2$ for which we obtain quasi-periodic trajectories with the frequencies $\omega_1, \omega_2$. Thus, in method (ⅰ) the series are parameterized by the amplitudes $A_1, A_2$, which can be selected at the beginning of the construction, while in method (ⅱ) the solutions are parameterized by the frequencies $\omega_1, \omega_2$, which are the parameters to select at the beginning of the construction. Also, in the latter case the series inverse to (1.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 (1.3). Thus, with all frequencies of the problem fixed in advance, establishing the convergence of the inverse series (1.5) suffices to answer the question posed at (ⅱ).

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 (ⅰ) we obtain divisors of the form $k_1\omega_{0, 1}+k_2\omega_{0, 2}$, with $(k_1, k_2)\in\mathbb{Z}^2$, $|k_1|+|k_2|\neq 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 (ⅰ) we introduce the parametrization (modulo two unimportant phases) $\varphi_1 = \omega_1 t$, $\varphi_2 = \omega_2 t$, and after introducing the series expressions (1.2) and (1.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 $\Theta_{1, i}(\varphi_1, \varphi_2)$, $\Theta_{2, i}(\varphi_1, \varphi_2)$ contain trigonometric terms in the angles $\varphi_1, \varphi_2$ (see ^[15], section 4).

- In scheme (ⅱ), instead, we obtain divisors of the form $k_1\omega_1+k_2\omega_2$, i.e., depending on the (fixed) pre-selected new frequencies $\omega_1, \omega_2$. This follows from the fact that the linear non-homogeneous equations to solve are now of the form (see ^[13]):

again with functions $\Phi_{1, i}(\varphi_1, \varphi_2)$, $\Phi_{2, i}(\varphi_1, \varphi_2)$ containing trigonometric terms in the angles $\varphi_1 = \omega_1 t, \varphi_2 = \omega_2 t$. Note that since the divisors depend on the new frequencies $\omega_1, \omega_2$, choosing non-resonant values for the latter permits the formal construction to proceed; this, even when the unperturbed frequencies $\omega_{0, 1}, \omega_{0, 2}$ are, instead, resonant.

As regards convergence, in the case (ⅰ) Poincaré (^[19], Ch.Ⅸ) already emphasizes that the Lindstedt series with divisors depending on the original harmonic frequencies $\omega_{0, 1}, \omega_{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 (ⅰ). 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 (ⅰ) 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 (ⅱ) is actually possible by diagrammatic methods via the following theorem ^[4]:

Theorem 1.1 (^[4]). Consider the $N$ coupled oscillator equations

where $\varepsilon$ is a real parameter, ${\mathit{\boldsymbol{f}}}({\mathit{\boldsymbol{x}}}, \varepsilon)$ is real analytic at ${\mathit{\boldsymbol{x}}} = 0\, , \, \varepsilon = 0$ and at least quadratic in ${\mathit{\boldsymbol{x}}}$, and the frequency vector ${\mathit{\boldsymbol{\omega}}}$ is diophantine. Let

be a solution of (1.8) for any choice of the complex constants $c_j$ and for $\varepsilon = 0$. Let $\Gamma({\mathit{\boldsymbol{c}}}) = \max(|c_1|, \ldots, |c_N|, 1)$. Then, there exists a positive constant $\eta_0$ and a function ${\mathit{\boldsymbol{\eta}}}(\varepsilon, {\mathit{\boldsymbol{c}}})$ holomorphic in the domain $|\varepsilon|\Gamma^3({\mathit{\boldsymbol{c}}})\leq\eta_0$, real for real $\varepsilon$, such that the system

admits a solution of the form

holomorphic in the domain $|\varepsilon|\Gamma^3({\mathit{\boldsymbol{c}}})e^{3|{\mathit{\boldsymbol{\omega}}}||\Im t|}\leq \eta_0$ and real for real $\varepsilon, t$. The constants ${\mathit{\boldsymbol{A}}}_{\mathit{\boldsymbol{\nu}}}$ are $\mathcal{O}(\varepsilon)$, except for the constants $A_{1, 0, ..., 0}$, $A_{0, 1, ..., 0}$, $A_{0, 0, ..., 1}$, which are equal to $c_1$, $c_2 , \ldots, 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:

where the frequencies $\Omega_k$, $k = 1, \ldots, N$ are given in terms of the FPU normal mode spectrum $\Omega_k = 2\sin(k\pi/(2(N+1))$, the function $F$ can be cubic or quartic in the variables $Q_k$, and the perturbation $\varepsilon$ 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(\omega_q t+\phi_q)$, and $Q_k(t) = 0$ for $k\neq q$. For the frequencies $\omega_q$ we obtain series expressions of the form

Then, for $\varepsilon$ sufficiently small, the Lindstedt method (ⅱ) above allows to represent the q-breathers via the Fourier expansion

where $\tilde{f}_{k, m} = {\mathcal{O}}(\varepsilon^{p(k, q, m)})$, with integer exponent $p(k, q, m)\geq 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 ${\mathcal{R}}_q = (q_1, \ldots q_M)\in \{1, 2, ..., N\}^M$. The corresponding Fourier representation of these special solutions can again been computed using Lindstedt series, and it obtains the form

with ${\mathit{\boldsymbol{m}}} \equiv (m_1, \ldots, m_M)\in\mathbb{Z}^M$, ${\mathit{\boldsymbol{A}}}\equiv(A_{q_1}, \ldots, A_{q_M})$, ${\mathit{\boldsymbol{\omega}}}\equiv(\omega_{q_1}, \ldots, \omega_{q_M})$, ${\mathit{\boldsymbol{\phi}}}\equiv(\phi_{q_1}, \ldots, \phi_{q_M})$, and $\tilde{f}_{k, {\mathit{\boldsymbol{m}}}} = {\mathcal{O}}(\varepsilon)$ for all $k = 1, \ldots 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 full-dimensional 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 ${\mathit{\boldsymbol{\omega}}}_0$. Some results and detailed proofs are deferred to the Appendix.

2.   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 = \sqrt{2J}\sin(q)\,$, $p = \sqrt{2 J}\cos(q)\,$, we obtain

Let $J_0$ be the label of a given torus (periodic orbit) of the harmonic oscillator model ${\mathcal{H}}_0$. Consider a real neighborhood ${\cal D}_\varepsilon = \{J = J_0+p\; \mbox{with}\; |p| < D_\varepsilon\}$, where $D_\varepsilon = {\cal O}(\varepsilon)$. We will illustrate four different perturbative methods to treat the dynamics in the phase-space neighborhood $\mathbb{T}\times{\cal D}_\varepsilon$: these are ⅰ) a Birkhoff normal form construction, with ⅱ) its analog in terms of Lindstedt series, ⅲ) a Lindstedt series exhibiting the torus-fixing property of the Kolmogorov method, and, finally ⅳ) the normal form analogue of ⅲ). 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 ${\mathcal{H}}_0$.

2.1.   Birkhoff normal form

Setting $J = J_0+p\,$ the Hamiltonian takes the form (apart from a constant)

A 'Birkhoff normal form' for the Hamiltonian (2.3) can be computed by introducing a canonical transformation eliminating the angle $q\,$. Writing the Hamiltonian (2.3) as

where $Z_0 = \omega_0 p$ and ${f_1 = \sum\limits_{k, l} c_{k, l}\, p^k e^{i l q}}\,$ we will define a Lie generating function $\chi^{(1)}(q, p)\,$ bringing this Hamiltonian to normal form up to terms ${\cal O}(\varepsilon)$. In the standard procedure, it is sufficient to set $\chi^{(1)} = X^{(1)}$, where $X^{(1)}$ solves the homological equation $L_{X^{(1)}}Z_0+\varepsilon f_1 = \varepsilon \zeta_1\,$, with $\zeta_1 = < f_1 > _q$ and $L_X$ denoting the Poisson bracket operator $L_X = \{\cdot, 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 $\chi^{(1)}$. Consider the canonical transformation $(q, p)\rightarrow (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 $\left(q(q^{(r)} = 0, p^{(r)} = 0), p(q^{(r)} = 0, p^{(r)} = 0)\right) = (0, 0)+{\cal O}(\varepsilon^{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 $\{\omega_0 p, K^{(n)}\, q\} = -\omega_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)

where

Omitting details, the formulas obtained after two normalization steps as above are the following: the Hamiltonian takes the form

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 $({\tilde{q}}, {\tilde{p}})$ in any function of the form $F^{(r)}({\tilde{q}}, {\tilde{p}})$ imply the new canonical variables computed after $r$ normalization steps. Then, up to order $2$ in $\varepsilon$ we obtain:

where the remainder $R^{(2)}$ is ${\cal O}(\varepsilon^3)$. The Hamiltonian $Z^{(2)}$ can now be used to analytically compute $\varepsilon^2-$precise solutions to the equations of motion in the variables $({\tilde{q}}, {\tilde{p}})$. The equations of motion are

fixing the initial condition ${\tilde{p}}(0) = {\tilde{q}}(0) = 0\,$ yields the solution

This can be back-transformed to the solution in the original variables. We have

Substituting the solutions $({\tilde{q}}(t), \, {\tilde{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:

(ⅰ) 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 $\omega_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 ${\mathit{\boldsymbol{m}}}\cdot{\mathit{\boldsymbol{\omega_0}}}$, where ${\mathit{\boldsymbol{m}}}\in\mathbb{Z}^N$, $|{\mathit{\boldsymbol{m}}}|\neq 0$, and ${\mathit{\boldsymbol{\omega_0}}}$ is the $N-$vector of the unperturbed frequencies of the problem. This implies that the method may formally proceed only when the vector ${\mathit{\boldsymbol{\omega_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).

(ⅱ) 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 ${{\mathit{\boldsymbol{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) $\omega$ 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 ${\cal O}(\varepsilon^2)$ in Eq (2.4).

(ⅲ) Any choice is possible to impose on the initial phase $q(0)$, altering slightly the definition of the constants $S^{(n)}$.

2.2.   Lindstedt solution analogous to Birkhoff

Analytical solutions in the form of formal series obtained like (2.6) 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. Ⅸ) 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 (2.6) 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 (2.2) are

We perform the following steps:

Step 1: time re-parametrization. Set $\varphi = \omega\, t$ where $\omega$ 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 $\omega$ as

Note that the corrections $a_i$ $i = 1, 2, \ldots$ 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, \ldots$ are considered functions of the phase $\varphi = \omega\, 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 $\varphi$ (or, more generally, of a set of angles ${\mathit{\boldsymbol{\varphi}}}\in\mathbb{T}^N$, with ${\mathit{\boldsymbol{\varphi}}} = {\mathit{\boldsymbol{\omega}}}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 ${\mathit{\boldsymbol{\omega_0}}}$. This is made clear in the next step.

Step 4: splitting of the equations of motion in powers of $\varepsilon$ and iterative solution. In our example, from the definition of $\varphi$ we have:

Substituting (2.8) and (2.9) in the equations of motion (2.7a) and (2.7b) leads to the following expressions (up to order $2$ in $\varepsilon\,$)

Now, in order to iteratively determine the functions $q_i(\varphi), J_i(\varphi)$, 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(\varphi) = \varphi$ and $J_0(\varphi) = 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

with $\varphi = \omega\, t$ and ${\omega = \omega_0+\frac{3 \, \varepsilon\, J_0}{4}-\frac{69 \, \varepsilon^2 \, J_0^2}{64 \, \omega_0}}\,$. Comparing the result with Eqs (2.4) and (2.6), 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 (2.10). After the expansion of the frequency, the differential equations to solve at subsequent steps all involve the operators $\omega_0d/d\varphi$. Hence, all divisors appearing in the series terms are in terms of the unperturbed frequency $\omega_0$.

2.3.   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.^* 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 ${\mathcal{H}}_0({\bf{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).

^*The examples treated in this section are obvious exceptions, since, we deal, for simplicity, with 1DOF systems exhibiting no small divisors. Small divisors appear, instead, for $N > 1$.

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 ${\omega}$. 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 ${\omega(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 (2.7a) and (2.7b) represented as a trigonometric series and evolving according to a given pre-selected frequency $\omega$. Inverting the expansion (2.8) 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 (2.7a) and (2.7b) the equations of motion lead now to the expressions (up to order $2$ in $\varepsilon\,$)

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:

which corresponds to the choice of the initial condition $q_0(0) = 0$ and $J_0(0) = J_0\,$. Then, at first order we have

implying $a_1 = -3 \, J_0/4\,$ and

Note that the integration constants in (2.16) 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 (2.17) above with those (Eq (2.12)) found in the case of the original Lindstedt method. We have the following remarks:

(ⅰ) From second order and beyond the coefficients in front of the harmonics $\cos(m\omega t), \sin(m\omega t)$ for the same $m$ are not all equal in the two solutions.

(ⅱ) Most importantly, the structure of the divisors in the two solutions is different. In the torus fixing case (Eq (2.17)), all divisors involve the corrected frequency, $\omega$, instead of the original frequency $\omega_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 ${\mathit{\boldsymbol{m}}}\cdot{\mathit{\boldsymbol{\omega}}}$. Thus, the method can proceed even when the original frequencies of the $N$ linear oscillators ${\mathit{\boldsymbol{\omega_0}}}$ are resonant, as long as we impose a non-resonant detuning in the adopted values for the final frequencies ${\mathit{\boldsymbol{\omega}}}$ on the torus.

(ⅲ) On the other hand, the value of the initial condition $J_0$ leading to motion on the torus with frequency $\omega$ 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.

(ⅳ) At the end of the process, however, $J_0$ can be estimated by reversing the series (2.13), 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

where the functions $a_{n, j}(J_{0, 1}, \ldots, J_{0, N})$ have all been specified iteratively up to a maximum truncation order in $n$, and they are polynomial in ${\mathit{\boldsymbol{J_0}}}\equiv(J_{0, 1}, \ldots, J_{0, N})$. Thus, the series (2.18) can be formally inverted, yielding

where the functions $P_{n, j}({\mathit{\boldsymbol{\omega}}}-{\mathit{\boldsymbol{\omega_0}}})$ are polynomial of degree $n$ in the differences $({\mathit{\boldsymbol{\omega}}}-{\mathit{\boldsymbol{\omega_0}}})$. Note that for the inverse series to converge we must require the difference $|({\mathit{\boldsymbol{\omega}}}-{\mathit{\boldsymbol{\omega_0}}})|$ to be smaller than $\varepsilon$, a fact which limits how far we can detune ${\mathit{\boldsymbol{\omega}}}$ from ${\mathit{\boldsymbol{\omega_0}}}$ in order to be able to specify the corresponding initial condition ${\mathit{\boldsymbol{J_0}}}$. An obvious choice is $|({\mathit{\boldsymbol{\omega}}}-{\mathit{\boldsymbol{\omega_0}}})| = \mathcal{O}(\varepsilon^2)$. At any rate, we emphasize that the convergence of the inverse series (2.18) is an open problem, crucial to the applications.^†

^†In the lack of a known answer to this problem, we can adopt two different attitudes: we may fix ${\mathit{\boldsymbol{J_0}}}$ and claim that our series answer the question of what was the value of ${\mathit{\boldsymbol{\omega_0}}}$ in the original Hamiltonian corresponding to the motion on the particular torus with frequency ${\mathit{\boldsymbol{\omega}}}$ (see ^[20]). On the other hand, in real world applications the only certainty we have is about our initial model (!), i.e., ${\mathit{\boldsymbol{\omega_0}}}$, so, we must devise a method to propose values of ${\mathit{\boldsymbol{\omega}}}$ for which the method will probably turn to converge. In the case of the FPU problem, a numerical procedure to move within the space of the parameters ${\mathit{\boldsymbol{\omega}}}$ so at to choose plausible values was proposed in ^[3].

2.4.   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 $\omega_0$ in the Hamiltonian (2.3) with the series (2.13). 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) ${\mathcal{H}}^{(2)}(q, p) = Z^{(2)}(q, p)+R^{(2)}(q, p)$ where the remainder $R^{(2)}$ is $\mathcal{O}(\varepsilon^3)$, while $Z^{(2)}(q, p)$ has the form $Z^{(2)}(q, p) = \omega\, p+ \varepsilon\, R_1(q, p)+\varepsilon^2\, R_2(q, p)$ with both $R_i(q, p) = {\mathcal{O}}(||\, {p} \, ||^2)$ $i = 1, \, 2$. To this end:

- First order: we fix $a_1$ so that the linear term $\varepsilon a_1 p$ acts as counterterm for the term $\varepsilon (3J_0/4)p$. This provides, in the twistless case, a process by which the frequency $\omega$ can be kept fixed (in the usual twist case, instead, this would have been accomplished by exploiting the Hessian matrix of ${\mathcal{H}}_0\,$). Formally, we require that

leading to $a_1 = -\frac{3\, J_0}{4}$. Now, we insert this expression for $a_1$ in the Hamiltonian (2.20), leading to

We can now eliminate the $O(\varepsilon)$ 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

leading to $L_{X^{(1)}(q)}\left(\omega\, p\right) + h_{1, 0}^{(0)} = \langle {h_{1, 0}^{(0)}} \rangle_{{q}}\,$, that is

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 $\chi_{1}^{(1)}(q) = X^{(1)}(q) + K^{(1)}q\,$, where $K^{(1)} = 3 \, J_0^2/(8\, \omega)$ is a constant. The general rules for the determination of the constants $K^{(j)}$ will be discussed in Section 3.

Using the generating function $\chi_{1}^{(1)}(q)$ we obtain the intermediate Hamiltonian

given by

We will now eliminate the trigonometric terms linear in the momentum in $\widehat{{\mathcal{H}}}^{(1)}$. These are

and can be eliminated by a generating function of the form $\chi_{2}^{(1)}(q, p) = \tilde{\chi}_{2}^{(1)}(q, p) + S^{(1)}\, p$ satisfying the homological equation $L_{\chi_{2}^{(1)}(q, p)}\left(\omega\, p\right)+\widehat{h}_{1, 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). Then

The new Hamiltonian is

and it is in Kolmogorov normal form up to order $\varepsilon\,$.

Repeating the same procedure at second order we arrive at the formulas

leading to the Hamiltonian ${\mathcal{H}}^{(2)}({\tilde{q}}, {\tilde{p}}) = Z^{(2)}({\tilde{q}}, {\tilde{p}})+R^{(2)}({\tilde{q}}, {\tilde{p}})$ with the Kolmogorov normal form part

($({\tilde{q}}, {\tilde{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({\tilde{q}}, {\tilde{p}}) = {\mathcal{O}}(||\, {{\tilde{p}}} \, ||^2)$ $j = 1, \, 2\,$, the equations of motions under the Hamiltonian $Z^{(2)}$ are

The torus ${\tilde{p}}(t) = 0\,$, ${\tilde{q}}(t) = \omega\, t$ (where we chose ${\tilde{q}}(0) = 0\,$) is a solution of this system. This can be back-transformed in the original variables using

Substituting the solution $({\tilde{q}}(t) = \omega\, t, \, {\tilde{p}}(t) = 0)$ in the previous expressions, we readily deduce the solution in the original variables:

Also, using the computed expressions of $a_1$ and $a_2$ (Eqs (2.21) and (2.24)), 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 (2.17)).

2.5.   Comparisons and numerical tests

In order to better visualize the differences between the methods (ⅰ) and (ⅱ) (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 $\varepsilon^2\,$.

From Table 1 above, we observe that the solutions $q(t)$ and $J(t)$ obtained by methods (ⅰ) and (ⅱ) have the same form up to order ${\mathcal{O}}(\varepsilon)\,$, except for the fact that in method (ⅰ) all divisors depend on the frequency $\omega_0$ rather than $\omega\,$. On the other hand, the coefficients found by the two methods start differing from the order $\varepsilon^2$ and beyond. Moreover, the series for the frequency $\omega$ 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\geq2\,$). To this end, we report below a comparison between the numerical solutions and the analytical ones, obtained by the methods (ⅰ) and (ⅱ) in the example of the Hamiltonian (2.1). To produce the Lindstedt series (ⅱ) in this example we work as follows: fixing the initial and final frequency $\omega_0$ and $\omega\,$, we reverse the series (2.8) according to

where $d\omega = \omega-\omega_0\,$ and $f^{-1}$ denotes the series inverse to $f$. Then, having specified $J_0$ through the inverse series (2.27), we compute all numerical coefficients in the Lindstedt series (ⅰ) and (ⅱ) up to order $\varepsilon^4\,$. We analyze three different cases:

$1)$ $\varepsilon = 1\,$, $\omega_0 = 1$ and $\omega = 1.002\,$;

$2)$ $\varepsilon = 1\,$, $\omega_0 = 1$ and $\omega = 1.02\,$;

$3)$ $\varepsilon = 1\,$, $\omega_0 = 1$ and $\omega = 1.2\,$.

Let us note that all results are rescalable to different choices of $\varepsilon\,$. Reversing the series (2.8) as prescribed by (2.27), we obtain, respectively, the following values for the amplitude (initial action datum): $J_0 = 0.0026769 \, ,$ $0.0277048\, , 0.383509\,$. Substituting these values of $\varepsilon\,$, $\omega\,$, $\omega_0$ and $J_0$ in the solutions $q(t)$ and $J(t)$ we obtain the solutions, as functions of $t\,$, produced by methods (ⅰ) and (ⅱ). At the same time, it is possible to integrate the equations of motion (2.7a) and (2.7b) to produce a numerical solution, starting from the initial conditions $q(0) = 0\,$, $J(0) = J_0\,$.

Figure 1 shows the difference (in $\mathrm{Log}_{10}\,$ scale) between the semi-analytical solutions obtained by methods (ⅰ) and (ⅱ) 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\omega$ is of order $10^{-3}$ and raising up to two orders of magnitude when $d\omega$ becomes of the order of unity.

2.6.   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\in\mathbb{R}$. This changes, however, in more general models in which the Hamiltonian is of the form

where the development in powers of the variables ${\mathit{\boldsymbol{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 $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 $\alpha-$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 ${\mathcal{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.

3.   KAM algorithm for isochronous systems

We now give the Kolmogorov algorithm for generic isochronous systems with Hamiltonian

where ${\mathit{\boldsymbol{q}}}\in\mathbb{T}^{n}\,$, ${\mathit{\boldsymbol{J}}}\in\mathcal{B}\subseteq\mathbb{R}^{n}\,$, $\tilde{h}_{i, \, s} = \mathcal{O}(||{\mathit{\boldsymbol{J}}}||^{\frac{s}{2}})\,$. We assume that the Hamiltonian (3.1) has a half-integer power dependence on ${\mathit{\boldsymbol{J}}}\,$, i.e., admits the expansion, for $n\rightarrow \infty$, of the truncated series

where $\tilde{h}^{(k)}$ are the $k-$th derivatives of $\tilde{h}$ with respect to ${\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{J}}}-{\mathit{\boldsymbol{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 ${\mathit{\boldsymbol{\omega}}}$ fixed in advance, given by

where the parameters ${\mathit{\boldsymbol{J_{0}}}}$, whose values are to be specified in the end of the process, give the initial conditions for ${\mathit{\boldsymbol{J}}}$ of a trajectory on the torus with frequencies ${\mathit{\boldsymbol{\omega}}}\,$. To this end, along with the normalizing canonical transformation, the algorith computes 'on the go' the functions ${\mathit{\boldsymbol{a}}}_{i}({\mathit{\boldsymbol{J}}}_{0})$ ('counter-terms').

Substituting the series (3.4) in the Hamiltonian (3.3) we arrive at:

In the following we use the notation

where $h_{i, 0}^{(k)}$ is independent from ${\mathit{\boldsymbol{p}}}\,$, $h_{i, 1}^{(k)}$ linear in ${\mathit{\boldsymbol{p}}}$ and the remaining sum is ${\mathcal{O}}(||\, {{\mathit{\boldsymbol{p}}}} \, ||^2)\,$. To facilitate reading, the indices $(i, j, k)$ used in all subsequent expressions refer to

We now have the following

Proposition 1. Assume the vector ${\mathit{\boldsymbol{\omega}}}$ is non-resonant. There exist Lie generating functions $\chi_1^{(r)}$, $\chi_2^{(r)}$ such that, after $r$ normalization steps, the Hamiltonian $(3.5)$ is given by the formal series:

where $h_{j}^{(r)} = \sum_{i\geq 0} h_{j, i}^{(r)}\, \forall\, j\geq r+1\,$ and $h_{j}^{(r)} = \sum_{i\geq 2} h_{j, i}^{(r)}\, \in\mathcal{O}(||{\mathit{\boldsymbol{p}}}||^{2})\, \forall\, j = 1, \ldots, r$, with

Given the Fourier series ${h_{r, 0}^{(r-1)}({\mathit{\boldsymbol{q}}}) = \sum\limits_{{\mathit{\boldsymbol{k}}}\in\mathbb{Z}^n} c_{{\mathit{\boldsymbol{k}}}}^{(r-1)}e^{i {\mathit{\boldsymbol{k}}}\cdot{\mathit{\boldsymbol{q}}}}}$ and ${h_{r, 1}^{(r-1)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}}) = \sum\limits_{|{\mathit{\boldsymbol{j}}}| = 1}\sum\limits_{{\mathit{\boldsymbol{k}}}\in\mathbb{Z}^n}c_{{\mathit{\boldsymbol{j}}}{\mathit{\boldsymbol{k}}}}^{(r)}{\mathit{\boldsymbol{p}}}^{{\mathit{\boldsymbol{j}}}}e^{i {\mathit{\boldsymbol{k}}}\cdot{\mathit{\boldsymbol{q}}}}}\,$, the generating functions $\chi_{1}^{(r)}$, $\chi_{2}^{(r)}$ are defined by

Proof of the proposition. The generic $r$-th iterative step of the algorithm is defined as follows: after $r-1$ steps, the hamiltonian (3.5) has the form:

where $C_{1}\, \ldots\, C_{r-1}$ are constants and $h_{i}^{(r-1)}\in {\mathcal{O}}(||\, {{\mathit{\boldsymbol{p}}}} \, ||^{2})\,$ $\forall\, i\, = \, 1, \ldots, r-1\,$. Taking into account the notation (3.7), we re-define the quantity $h_{i}^{(k)}$ according to

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 $\chi_{1}^{(r)}$ is defined as $\chi_{1}^{(r)}({\mathit{\boldsymbol{q}}}) = X^{(r)}({\mathit{\boldsymbol{q}}})+{\mathit{\boldsymbol{K}}}^{(r)}\cdot{\mathit{\boldsymbol{q}}}\,$, with ${\mathit{\boldsymbol{K}}}^{(r)}$ an appropriate constant vector defined below, and $X^{(r)}$ is defined through the homological equation:

where $\langle {.} \rangle$ denotes the mean over ${\mathit{\boldsymbol{q}}}\,$. The function $X^{(r)}({\mathit{\boldsymbol{q}}})$ eliminates all terms of order ${\mathcal{O}}(\varepsilon^r)\,$ depending only on the angles ${\mathit{\boldsymbol{q}}}\,$ in the Hamiltonian ${\mathcal{H}}^{(r-1)}$. Writing $h_{r, 0}^{(r-1)}({\mathit{\boldsymbol{q}}})$ in the Fourier form ${h_{r, 0}^{(r-1)}({\mathit{\boldsymbol{q}}}) = \sum\limits_{{\mathit{\boldsymbol{k}}}\in{\mathbb{Z}}^n} c_{{\mathit{\boldsymbol{k}}}}^{(r-1)}e^{i {\mathit{\boldsymbol{k}}}\cdot{\mathit{\boldsymbol{q}}}}}$ from (3.11) we find:

We then impose the condition that the terms of order ${\mathcal{O}}(\varepsilon^{r})$ linear in ${\mathit{\boldsymbol{p}}}$ have zero

This specifies ${\mathit{\boldsymbol{a}}}_{r}$ as a function of ${\mathit{\boldsymbol{J}}}_{0}\,$ via Eq (3.9).

Finally, we specify the constant vector ${\mathit{\boldsymbol{K}}}^{(r)}\,$. To this end, we impose the condition that the solution in ${\mathit{\boldsymbol{p}}}$ has the form $\sum_{{\mathit{\boldsymbol{k}}}} A_{{\mathit{\boldsymbol{k}}}}\left(\cos({\mathit{\boldsymbol{k}}}\cdot{\mathit{\boldsymbol{q}}})-1\right)\,$, so that, at time $t = 0$, we have $p(0) = 0$. Writing (3.12) as

the poisson bracket of $\chi_{1}^{(r)}({\mathit{\boldsymbol{q}}})$ with the i-th component $p_i\,$ ($1 \leq \, i \, \leq n\,$) of the vector ${\mathit{\boldsymbol{p}}}$ yields the following expressions for the i-th components $K_i^{(r)}$ and $k_i$ of the vectors ${\mathit{\boldsymbol{K}}}^{(r)}$ and ${\mathit{\boldsymbol{k}}}\,$:

Hence

i.e.,

Finally, we compute $\widehat{{\mathcal{H}}}^{(r)}$ in (3.10) as:

where $\widehat{h}_{j}^{(r)} = \sum_{i\geq 0} \widehat{h}_{j, i}^{(r)}\, \forall\, j\geq r\,$ and $\widehat{h}_{j}^{(r)} = \sum_{i\geq 2} \widehat{h}_{j, i}^{(r)}\in{\mathcal{O}}(||\, {{\mathit{\boldsymbol{p}}}} \, ||^{2})\,$ $\forall\, j = 1, \ldots, r-1\,$, where:

● Second half step (of r-th): We compute

with a generating function $\chi_{2}^{(r)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}})$ linear in ${\mathit{\boldsymbol{p}}}$. Setting $\chi_{2}^{(r)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}}) = \tilde{\chi}_{2}^{(r)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}})+ {\mathit{\boldsymbol{S}}}^{(r)}\cdot{\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Y}}}^{(r)}({\mathit{\boldsymbol{q}}})\cdot{\mathit{\boldsymbol{p}}}+ {\mathit{\boldsymbol{S}}}^{(r)}\cdot{\mathit{\boldsymbol{p}}}\,$. The function $\tilde{\chi}_{2}^{(r)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}})$ satisfies the homological equation

Setting ${\widehat{h}_{r, 1}^{(r)}({\mathit{\boldsymbol{q}}}, {\mathit{\boldsymbol{p}}}) = \sum\limits_{|{\mathit{\boldsymbol{j}}}| = 1}\sum\limits_{{\mathit{\boldsymbol{k}}}\in{\mathbb{Z}}^n}c_{{\mathit{\boldsymbol{j}}}{\mathit{\boldsymbol{k}}}}^{(r)}{\mathit{\boldsymbol{p}}}^{{\mathit{\boldsymbol{j}}}}e^{i {\mathit{\boldsymbol{k}}}\cdot{\mathit{\boldsymbol{q}}}}}$ the solution of (3.27) is

We finally compute the constant vector ${\mathit{\boldsymbol{S}}}^{(r)}\,$ by the condition ${\mathit{\boldsymbol{q}}}(0) = 0$ at the time $t = 0$. By the poisson bracket

we obtain

The Hamiltonian ${\mathcal{H}}^{(r)}$ (Eq (3.26)) is:

where $h_{j}^{(r)} = \sum_{i\geq 0} h_{j, i}^{(r)}\, \forall\, j\geq r+1\,$ and $h_{j}^{(r)} = \sum_{i\geq 2} h_{j, i}^{(r)}\, \in{\mathcal{O}}(||\, {{\mathit{\boldsymbol{p}}}} \, ||^{2})\, \forall\, j = 1, \ldots, r\,$, where:

Equation (3.36), using homological equation (3.27), can be written equivalently as:

where, from (3.36) we have

Second part of the proof: using the formulas (3.19)–(3.25) and (3.32)–(3.38), we can express each term in the function $h^{(r)}$ in terms of the functions $h^{(r-1)}\,$ instead of $\widehat{h}^{(r)}\,$. From Eq (3.32) we have

where we used Eqs (3.19), (3.20) and (3.21), respectively. Using (3.25) we have

Equation (3.38) can now be written in the form

$\forall \: k\geq r+1\, , \: i\geq 2\,$. The indices $i, j, k$ in the first term of the above equation satisfy the relation

where we have used the inequalities $j\leq \lfloor (k-1)/r \rfloor -1$ and $\lfloor k-1 \rfloor \leq k-1\,$. Then Eq (3.45) ensures that the first term of (3.44) satisfies the definition (3.24). For the second term in (3.44), we have the following useful (also in the sequel) remark: we can write $k = mr+f\,$, where $0\leq f\leq r-1\,$ and $m\in\mathbb{N}\,$. Thus

Then

and

We can conclude that $1\leq k-\lfloor (k-1)/r \rfloor r \leq r\,$, implying that the definition (3.25) holds for the second term in (3.44). We can then write (3.44) as^‡

^‡Observe that, due to (3.47) and (3.48), we have that if $j = \lfloor (k-1)/r \rfloor\,$, then $0\leq \lfloor(k-jr-1)/r \rfloor = \lfloor(k-\lfloor(k-1)/r \rfloor r -1)/r \rfloor \leq\lfloor (r-1)/r \rfloor = 0\,$, i.e., $s = 0\,$. This allows to join all terms in a single sum.

$\forall \: k\geq r+1\, , \: i\geq 2\,$.

Equation (3.33) can be analyzed similarly as above, by splitting it in three different parts

$\forall \: k\geq 2r+1\,$. We study separately the relations satisfied by the indices $(j, k)$ of each of the three terms in the previous equality. Following (3.45), we have that

(since $j\leq \lfloor (k-1)/r \rfloor -2 \,$). Then, the definition (3.21) 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$.

(ⅰ) First case: $k$ is not a multiple of $r\,$, i.e.,

From (3.48) we have that $1\leq k-\lfloor \frac{k-1}{r} \rfloor r \leq r-1$ and, consequently, $r+1\leq k-\lfloor \frac{k-1}{r} \rfloor r + r\leq 2r-1\,$. Then, the definitions (3.19) and (3.20) hold, respectively, for the third and second term of (3.50). Thus, we can write Eq (3.50) as

$\forall \: k\geq 2r+1\,$, $k\neq mr\,$, $m\in\mathbb{N}\,$.

(ⅱ) Second case: $k$ is a multiple of $r\,$, i.e.,

From (3.47) we now have that $k-\lfloor \frac{k-1}{r} \rfloor r = r$ and, consequently, $k-\lfloor \frac{k-1}{r} \rfloor r +r = 2r\,$. Then, the definitions (3.19) and (3.21) hold, respectively, for the third and the second term of (3.50). Thus, we can write Eq (3.50) as

$\forall \: k\geq 2r+1\,$, $k = mr\,$, $m\in\mathbb{N}\,$. By the same argument, we can write Eq (3.35) as:

$\forall \: k\geq r+1\, , k\neq m\, r\, (m\in\mathbb{N})\,$. In view of the inequalities (3.45) and (3.48) we then readily find that the definitions (3.24) and (3.22) hold, respectively, for the first and second part of the previous equation, i.e., (3.52) leads to

$\forall \: k\geq r+1\, , k\neq m\, r\, (m\in\mathbb{N})\,$. Finally, recalling again (3.45) and (3.47), we can write Eq (3.36) as^§

^§In the last equality, due to (3.47), we have that if $j = \lfloor (k-1)/r \rfloor\,$, then $\lfloor(k-j r -1)/r \rfloor = \lfloor(k-\lfloor(k-1)/r \rfloor r -1)/r \rfloor = \lfloor (r-1)/r \rfloor = 0\,$, i.e., $s = 0\,$. Thus, also here we can join the terms in a single sum.

$\forall \: k\geq r+1\, , k = m\, r\, \, (m\in\mathbb{N})\,$ (where we used the definitions (3.24) and (3.23) for the first and second part of the sum). As before, we can then write the previous equation as

Use of Eq (3.23) and the homological equation (3.27) then concludes the proof.

4.   Concluding remarks

The focus of the present paper is twofold.

ⅰ) 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.

ⅱ) 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]).

Conflict of interest

The authors declare no conflict of interest.

A.   Example – 1DOF Hamiltonian with an odd (cubic) degree dependence on $x$

Consider the following one-degree of freedom Hamiltonian, with odd power dependence on the variable $x$

We pass to action-angle variables $(J, \, q)$ through the transformation $x = \sqrt{2J}\sin(q)\,$, $p = \sqrt{2 J}\cos(q)\,$, obtaining

Finally, we define the translation $J = J_0+p\,$ leading to (apart from constants)

Since the variable $J = J_0+p$ appears in the above Hamiltonian in half-integer powers, we expand the Hamiltonian in powers of the variable $p\,$ up to the same order as the maximum normalization order in $\varepsilon$. This proves to be sufficient since higher powers of $p$ only influence the process at powers (in the book keeping order) higher than the maximum normalization order. In particular, the following lemma can be easily proved:

Lemma 1. Let${\mathcal{H}}(q, J) = \omega_0\, J + \varepsilon\, h(q, J)\,$, where $h = {\mathcal{O}}(J^\frac{k}{2})\,$ $k\geq 3\,$. Then, for all $n\geq 1$ (with $n =$ number of normalization steps) $\chi_1^{(n)} = {\mathcal{O}}\left(J_0^{\frac{(k-2)n+2}{2}}\right)$ and $\chi_2^{(n)} = {\mathcal{O}}\left(J_0^{\frac{(k-2)n}{2}}\right)\,$.

Proof. For the proof, see the appendix B.

We will now illustrate the method by computing the direct (Lindstedt) and indirect (Kolmogorov) series up to order $3$ in $\varepsilon\,$ in the example above (the reason for reaching order $3$ instead of $2$ will become clear below). Starting from the Hamiltonian (A.3) we perform an expansion in $p$ leading to (apart from constants)

A.1.   Lindstedt series analogous to Kolmogorov

Following the same procedure as in subsection 2.3, we start with the equations of motion under the Hamiltonian (A.2)

Replacing, as before, the expressions

into the equations of motion and performing an expansion up to order $3$ in $\varepsilon\,$ (having fixed $\omega$) we compare terms of like orders in $\varepsilon$ in the l.h.s and r.h.s of (A.5a) and (A.5b). At order zero we have

where we fix the initial phase $q_0(0) = 0$ and $J_0(0) = J_0\,$. At order one, we find

Since no constant terms arise in $\dot{q_1}(t)$, we have that $a_1 = 0$. Then

yielding the constants $q_1(0) = J_1(0) = 0\,$. At order two we now get

To compensate for the constant term in $\dot{q_2}(t)\,$ we now set $a_2 = 5\, J_0/(6\, \omega)$. Then,

In a similar way, at order $3$ we find $a_3 = 0$, and the solutions

A.2.   Kolmogorov normal form

Starting from the Hamiltonian (A.4), we perform canonical transformations in order to bring the Hamiltonian into Kolmogorov normal form, i.e., ${\mathcal{H}}(q, p) = \omega\, p+\varepsilon\, R_1(q, \, p)+\varepsilon^2\, R_2(q, \, p)+\varepsilon^3\, R_3(q, \, p)$ where $R_{i}(q, \, p) = {\mathcal{O}}(||\, {p} \, ||^2)$ $i = 1, \, 2, \, 3$. Substituting, as in Section 2, the expression $\omega_0 = \omega + \varepsilon\, a_1+\varepsilon^2\, a_2+\varepsilon^3\, a_3\,$ in the Hamiltonian (A.4), we have

At first order, we have

implying that $a_1 = 0$ for the corresponding counterterm in the Hamiltonian (A.11). In order to eliminate the terms constant in the actions (depending only in the angle $q\,$) $\varepsilon\, h_{1, 0}$, given by

we define the generating function $X^{(1)}(q)$ such that ${L_{X^{(1)}(q)}(\omega\, p)+h_{1, 0}^{(0)} = \langle {h_{1, 0}^{(0)}} \rangle_{{q}}}\,$. Hence

In order to fix the initial value of $p$ at zero, we then set (as in Section 2) $\chi_{1}^{(1)}(q) = X^{(1)}(q) + K^{(1)}\, q\,$; in this case, the constant $K^{(1)} = 0\,$. We can now determine the intermediate Hamiltonian

Denoting by $\varepsilon\, \widehat{h}_{1, 1}^{(1)}$ the terms of $\widehat{{\mathcal{H}}}^{(1)}$ of order one in $\varepsilon$ and linear in $p\,$, we have

These terms are eliminated by the generating function $\chi_{2}^{(1)}(q, \, p) = \tilde{\chi}_{2}^{(1)}(q, \, p)+S^{(1)}\, p$ satisfying the homological equation $L_{\chi_{2}^{(1)}(q, \, p)}(\omega\, p)+\widehat{h}_{1, 1}^{(1)} = 0\,$. In order to fix the initial value of $q$ at zero, we readily find

and

Thus, the Hamiltonian

is now in Kolmogorov normal form up to first order in $\varepsilon\,$ (it is easy to check that the transformed Hamiltonian contains the normal form terms $- {\frac{3\, p^2 }{8\, \sqrt{2} \, \sqrt{J_0}}\, \sin (q)}+ {\frac{p^2}{8 \, \sqrt{2} \, \sqrt{J_0}}\, \sin (3\, q)}+ {\frac{p^3 }{16\, \sqrt{2}\, J_0^{3/2}}\, \sin (q)}- {\frac{p^3}{48\, \sqrt{2} \, J_0^{3/2}}\, \sin (3\, q)}\,$).

In a similar way we can proceed at orders 2 and 3, obtaining the formulas:

and

for the generating functions. The final Hamiltonian is

where $({\tilde{q}}, \, {\tilde{p}})$ indicate the new variables, and $R_j({\tilde{q}}, {\tilde{p}}) = {\mathcal{O}}(||\, {{\tilde{p}}} \, ||^2)$ $j = 1, 2, 3\,$. Thus, the torus ${\tilde{p}}(t) = 0, \, {\tilde{q}}(t) = \omega\, t\,$ is a solution for the equations of motion of this Hamiltonian. Using the Lie transformations, the solution in the original variables reads

Recalling also the computed values of $a_1\,$, $a_2\,$ and $a_3$ (Eqs (A.12)–(A.14)) we have also

Thus, we obtain the same solutions as by the Lindstedt method (Eq (A.10)).

We remark that also in this case the solutions produced by the Birkhoff normal form are equal to those produced by the Lindstedt method in the version 'analogous to Birkhoff'. However, if we compare these solutions with those produced by the KAM algorithm, we also note many differences from order $3$ and beyond. For completeness, we report in the following the solutions with the Birkhoff method:

B.   Proof of Lemma 1

Consider the Hamiltonian:

where $h = {\mathcal{O}}(J^{\frac{k}{2}})\,$ $k\in\mathbb{N}$, $k\geq 3\,$. Introducing the translation $J = J_0+p\,$, $h = {\mathcal{O}}((J_0+p)^{\frac{k}{2}})\,$, if $k$ is even we obtain a finite expression in the powers of $p$. If $k$ is odd, define $f(p) = \varepsilon\, (J_0+p)^{\frac{k}{2}}\,$ and suppose we want to compute the normalization steps by the Kolmogorov algorithm up to order $\varepsilon^n\,$ $n\geq 1\,$. The expansion in $p = 0\,$ up to order $n$ yields

Moreover, observe that $\chi_{1}^{(1)} \sim \varepsilon\, J_{0}^{\frac{k}{2}}$ and $\chi_{2}^{(1)} \sim \varepsilon\, J_0^{\frac{k}{2}-1}\,$. Now, we identify those terms whose Lie derivatives could modify the generating functions along the normalization process. We have $L_{\varepsilon\chi_{1}^{(1)}}^{n-1}\left(\varepsilon \, J_0^{\frac{k}{2}-n}p^n\right)$. Thus the term of higher degree in $p$ in the expansion of $f(p)\,$, denoted by $\simeq f^{(n)}(p)\,$, influences $h_{n, 1}^{(1)}$. In fact, since

after $n$ step a contribution stemming from the term $f^{(n)}(p)$ will appear in the generating function $\chi_{2}^{(n)}\,$. On the contrary, the term $L_{\varepsilon\chi_{1}^{(1)}}^{n-1}\left(\varepsilon \, f^{(n+1)}(p)\right)$ contributes to $h_{n, 2}^{(1)}$ (which is quadratic in the action $p$), while the term $L_{\varepsilon\chi_{1}^{(1)}}^{n}\left(\varepsilon \, f^{(n+1)}(p)\right)$ contributes to $h_{n+1, 1}^{(1)}$, which is of order $\varepsilon^{n+1}$. From (9.1) it follows that $h_{n, 1}^{(n)} = {\mathcal{O}}\left(J_0^{(k-2)n/2}\right)\,$. Hence, the dependence of $\chi_2^{(n)}$ on the parameter $J_0\,$ is in the power $\chi_2^{(n)}\sim\, J_0^{\frac{(k-2)n}{2}}\,$. Also

thus, again, $h_{n, 1}^{(n)} = {\mathcal{O}}\left(J_0^{(k-2)n/2}\right)\,$. Finally

implying $\chi_1^{(n)}\sim\, J_0^{\frac{(k-2)n+2}{2}}\,$. This concludes the proof of the lemma.