Kolmogorov algorithm
for isochronous Hamiltonian systems
R<EMAIL_ADDRESS>and C. Efthymiopoulos
Dipartimento di Matematica Tullio Levi-Civita, Università degli Studi di
Padova
via Trieste 63, 35121 Padova

###### Abstract

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.

##  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

$\mathcal{H}={1\over 2}(p_{x}^{2}+p_{y}^{2})+{1\over
2}\omega_{0,1}^{2}x^{2}+{1\over 2}\omega_{0,2}^{2}y^{2}+\varepsilon
P_{3}(x,y)$ (1)

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) under the form:

$x(t)=x_{0}(t)+\varepsilon
x_{1}(t)+\varepsilon^{2}x_{2}(t)+\ldots,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ y(t)=y_{0}(t)+\varepsilon
y_{1}(t)+\varepsilon^{2}y_{2}(t)+\ldots,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (2)

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:

$\displaystyle\omega_{1}(A_{1},A_{2})$ $\displaystyle=$
$\displaystyle\omega_{0,1}+\varepsilon\omega_{1,1}(A_{1},A_{2})+\varepsilon^{2}\omega_{2,1}(A_{1},A_{2})+\ldots,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\omega_{2}(A_{1},A_{2})$ $\displaystyle=$
$\displaystyle\omega_{0,2}+\varepsilon\omega_{1,2}(A_{2},A_{2})+\varepsilon^{2}\omega_{2,2}(A_{1},A_{2})+\ldots.$
(3)

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

$x_{0}(t)=A_{1}\cos(\omega_{1}t+\phi_{x0}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
y_{0}(t)=A_{2}\cos(\omega_{2}t+\phi_{y0})\leavevmode\nobreak\
\leavevmode\nobreak\ $ (4)

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

(ii) 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) 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 $\omega_{1},\omega_{2}$, invert the series (1) 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 (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 $\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) turn out to have the form (in the
cubic case)

$\varepsilon^{2}A_{1}^{2}=\sum_{i=1}^{\infty}C_{i,1}(\omega_{1}-\omega_{0,1})^{i},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\varepsilon^{2}A_{2}^{2}=\sum_{i=1}^{\infty}C_{i,2}(\omega_{2}-\omega_{0,2})^{i},$
(5)

for some constant coefficients $C_{i,1},C_{i,2}$ computable from the series
(1). 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}\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 (i) we introduce the
parametrization (modulo two unimportant phases) $\varphi_{1}=\omega_{1}t$,
$\varphi_{2}=\omega_{2}t$, and after introducing the series expressions (2)
and (1) to the equations of motion and separate terms of like orders we arrive
at equations (to be solved iteratively) of the form:

$\displaystyle\left(\omega_{0,1}{\partial\over\partial\varphi_{1}}+\omega_{0,2}{\partial\over\partial\varphi_{2}}\right)^{2}x_{i}+x_{i}=\Theta_{1,i}(\varphi_{1},\varphi_{2})$
$\displaystyle\left(\omega_{0,1}{\partial\over\partial\varphi_{1}}+\omega_{0,2}{\partial\over\partial\varphi_{2}}\right)^{2}y_{i}+y_{i}=\Theta_{2,i}(\varphi_{1},\varphi_{2})\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ i=1,2,\ldots$

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 (ii), 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]):

$\displaystyle\ddot{x}_{i}+\omega_{1}^{2}x_{i}=\Phi_{1,i}(\varphi_{1},\varphi_{2})$
$\displaystyle\ddot{y}_{i}+\omega_{2}^{2}y_{i}=\Phi_{2,i}(\varphi_{1},\varphi_{2})\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ i=1,2,\ldots$

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 (i) Poincaré ([19], Ch.IX) 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 (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]:

######  Theorem 1.1 ([4]).

Consider the $N$ coupled oscillator equations

$\ddot{x}_{j}+\omega^{2}_{j}x_{j}+f_{j}(x_{1},\ldots,x_{N};\varepsilon)+(\omega_{0,j}^{2}-\omega_{j}^{2})x_{j}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ j=1,\ldots,N$
(8)

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

$x^{(0)}_{j}(t)=c_{j}e^{i\omega_{j}t}+c_{j}^{*}e^{-i\omega_{j}t},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ j=1,\ldots,{N}$ (9)

be a solution of (8) for any choice of the complex constants $c_{j}$ and for
$\varepsilon=0$. Let $\Gamma({\bm{c}})=\max(|c_{1}|,\ldots,|c_{N}|,1)$. Then,
there exists a positive constant $\eta_{0}$ and a function
${\bm{\eta}}(\varepsilon,{\bm{c}})$ holomorphic in the domain
$|\varepsilon|\Gamma^{3}({\bm{c}})\leq\eta_{0}$, real for real $\varepsilon$,
such that the system

$\ddot{x}_{j}+\omega^{2}_{j}x_{j}+f_{j}(x_{1},\ldots,x_{{N}},\varepsilon)+\eta_{j}(\varepsilon,{\bm{c}})x_{j}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ j=1,\ldots,{N}$
(10)

admits a solution of the form

${\bm{x}}(t,\varepsilon,{\bm{c}})=\sum_{{\bm{\nu}}\in\mathbb{Z}^{{N}}}{\bm{A}}_{\bm{\nu}}\exp(i{\bm{\nu}}\cdot({\bm{\omega}}t))$
(11)

holomorphic in the domain
$|\varepsilon|\Gamma^{3}({\bm{c}})e^{3|{\bm{\omega}}||\Im t|}\leq\eta_{0}$ and
real for real $\varepsilon,t$. The constants ${\bm{A}}_{\bm{\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}$, …,
$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:

$\ddot{Q_{k}}+\Omega^{2}_{k}Q=\varepsilon F_{k}(Q_{1},\ldots,Q_{N})$ (12)

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

$\omega_{q}=\Omega_{q}+\Delta\omega_{q}(A_{q};\varepsilon),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \Delta\omega_{q}=\mathcal{O}(\varepsilon).$ (13)

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

$\displaystyle Q_{q}(t)$ $\displaystyle=$ $\displaystyle
A_{q}\cos(\omega_{q}t+\phi_{q})+\sum_{m=0}^{\infty}\tilde{f}_{q,m}(A_{q};\varepsilon)\cos[m(\omega_{q}t+\phi_{q})]$
(14) $\displaystyle Q_{k}(t)$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{\infty}\tilde{f}_{k,m}(A_{q};\varepsilon)\cos[m(\omega_{q}t+\phi_{q})],\leavevmode\nobreak\
\leavevmode\nobreak\ k\neq q,$

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

$\omega_{q_{i}}=\Omega_{q_{i}}+\Delta\omega_{q_{i}}(A_{q_{1}},A_{q_{2}},\ldots,A_{q_{M}};\varepsilon),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \Delta\omega_{q_{{i}}}=\mathcal{O}(\varepsilon),$ (15)

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

$\displaystyle Q_{k}(t)$ $\displaystyle=$ $\displaystyle
A_{k}\cos(\omega_{k}t+\phi_{k})+\sum_{|\mathbf{m}|=0}^{\infty}\tilde{f}_{k,{{\bm{m}}}}(\mathbf{A};\varepsilon)\cos[{\bm{m}}\cdot({\bm{\omega}}t+{\bm{\phi}})],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ k\in\mathcal{R}_{q}$ (16)
$\displaystyle Q_{k}(t)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\sum_{|\mathbf{m}|=0}^{\infty}\tilde{f}_{k,{{\bm{m}}}}(\mathbf{A};\varepsilon)\cos[{\bm{m}}\cdot({\bm{\omega}}t+{\bm{\phi}})],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ k\notin\mathcal{R}_{q}$

with ${\bm{m}}\equiv(m_{1},\ldots,m_{M})\in\mathbb{Z}^{M}$,
${\bm{A}}\equiv(A_{q_{1}},\ldots,A_{q_{M}})$,
${\bm{\omega}}\equiv(\omega_{q_{1}},\ldots,\omega_{q_{M}})$,
${\bm{\phi}}\equiv(\phi_{q_{1}},\ldots,\phi_{q_{M}})$, and
$\tilde{f}_{k,{\bm{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 ${\bm{\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

$\mathcal{H}(x,p)=\mathcal{H}_{0}+\varepsilon\mathcal{H}_{1}=\frac{\omega_{0}}{2}\left(p^{2}+x^{2}\right)+\varepsilon\frac{x^{4}}{4}\,.$
(17)

Using the harmonic oscillator action-angle variables $(J,\,q)$ with
$x=\sqrt{2J}\sin(q)\,$, $p=\sqrt{2J}\cos(q)\,$, we obtain

$\mathcal{H}(q,J)=\omega_{0}\,J+\frac{3\,\varepsilon}{8}J^{2}-\frac{\varepsilon}{2}J^{2}\cos(2q)+\frac{\varepsilon}{8}J^{2}\cos(4q)\,.$
(18)

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\leavevmode\nobreak\
\mbox{with}\leavevmode\nobreak\ |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 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 $\mathcal{H}_{0}$.

###  2.1 Birkhoff normal form

Setting $J=J_{0}+p\,$ the Hamiltonian takes the form (apart from a constant)

$\begin{split}\mathcal{H}(q,p)&=\omega_{0}\,p+\frac{3\,\varepsilon\,J_{0}^{2}}{8}+\frac{3\,J_{0}\,\varepsilon}{4}p+\frac{3\,\varepsilon}{8}p^{2}-\frac{\varepsilon\,J_{0}^{2}}{2}\cos(2q)-\varepsilon
J_{0}p\cos(2q)\\\
&\phantom{=}-\frac{\varepsilon}{2}p^{2}\cos(2q)+\frac{\varepsilon
J_{0}^{2}}{8}\cos(4q)+\frac{\varepsilon\,J_{0}}{4}p\cos(4q)+\frac{\varepsilon}{8}p^{2}\cos(4q)\,.\end{split}$
(19)

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

$\mathcal{H}(q,p)=Z_{0}+\varepsilon f_{1}\,,$

where $Z_{0}=\omega_{0}p$ and
$\displaystyle{f_{1}=\sum_{k,l}c_{k,l}\,p^{k}e^{ilq}}\,$ 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

$\chi^{(n)}(q,p)=X^{(n)}(q,p)+K^{(n)}\,q+S^{(n)}\,p\,,\quad n=1,\ldots,r$

where $K^{(n)}$, $S^{(n)}$ are constants possible to compute at every step by
requiring that

$\left(q(q^{(n)}=0,p^{(n)}=0),p(q^{(n)}=0,p^{(n)}=0)\right)=(0,0)+{\cal
O}(\varepsilon^{n+1}).$

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)

$\chi^{(1)}(q,p)=X^{(1)}(q,p)+K^{(1)}\,q+S^{(1)}\,p\,,$

where

$\displaystyle X^{(1)}(q,p)=\sum_{\begin{subarray}{c}l\neq 0\\\
k\end{subarray}}\frac{c_{k,l}}{i\,l\,\omega_{0}}\,p^{k}e^{ilq}\,,$
$\displaystyle K^{(1)}=-\sum_{l\neq 0}\frac{c_{0,l}}{\omega_{0}}\,,$
$\displaystyle S^{(1)}=-\sum_{l\neq 0}\frac{c_{1,l}}{i\,l\,\omega_{0}}\,.$

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

$\mathcal{H}^{(2)}(q^{(2)},p^{(2)})=\mathrm{exp}L_{\varepsilon^{2}\chi^{(2)}}\mathrm{exp}L_{\varepsilon\chi^{(1)}}\mathcal{H}\Big{|}_{\begin{subarray}{c}q=q^{(2)}\\\
p=p^{(2)}\end{subarray}}=\left(Z_{0}+\varepsilon
Z_{1}+\varepsilon^{2}Z_{2}+\sum_{i\geq
3}\varepsilon^{i}f_{i}\right)\\!\bigg{|}_{\begin{subarray}{c}q=q^{(2)}\\\
p=p^{(2)}\end{subarray}}\,.$

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:

$\displaystyle\mathcal{H}^{(2)}(\tilde{q},\tilde{p})=Z^{(2)}+R^{(2)}=\omega_{0}\,\tilde{p}+\frac{3\,\varepsilon\,J_{0}}{4}\,\tilde{p}+\frac{3\,\varepsilon\,{\tilde{p}}^{2}}{8}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega_{0}}\,\tilde{p}-\frac{51\,\varepsilon^{2}\,J_{0}}{64\,\omega_{0}}\,{\tilde{p}}^{2}-\frac{17\,\varepsilon^{2}}{64\,\omega_{0}}\,{\tilde{p}}^{3}+R^{(2)}\,$

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

$\displaystyle\begin{cases}&\displaystyle{\dot{\tilde{q}}=\frac{\partial\mathcal{H}^{(2)}}{\partial\tilde{p}}=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}+\frac{3\,\varepsilon}{4}\tilde{p}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega_{0}}-\frac{51\,\varepsilon^{2}\,J_{0}}{32\,\omega_{0}}\,\tilde{p}-\frac{51\,\varepsilon^{2}}{64\,\omega_{0}}\,\tilde{p}^{2}}\\\\[7.3194pt]
&\displaystyle{\dot{\tilde{p}}=-\frac{\partial\mathcal{H}^{(2)}}{\partial\tilde{q}}=0}\end{cases}\,;$

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

$\displaystyle\tilde{p}(t)=0\,,$
$\displaystyle\tilde{q}(t)=\omega\,t:=\varphi_{B}\,,$
$\displaystyle\omega=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega_{0}}\,.$
(20)

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

$\displaystyle
q=\mathrm{exp}L_{\varepsilon^{2}\,\chi^{(2)}}\mathrm{exp}L_{\varepsilon\,\chi^{(1)}}\,\tilde{q}\,,$
$\displaystyle
p=\mathrm{exp}L_{\varepsilon^{2}\,\chi^{(2)}}\mathrm{exp}L_{\varepsilon\,\chi^{(1)}}\,\tilde{p}\,.$
(21)

Substituting the solutions $(\tilde{q}(t),\,\tilde{p}(t))$ in the previous
expression, we find

$\begin{split}q(t)&=\varphi_{B}-\frac{\varepsilon\,J_{0}}{2\,\omega_{0}}\,\sin(2\varphi_{B})+\frac{31\,\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(2\,\varphi_{B})+\frac{\varepsilon\,J_{0}}{16\,\omega_{0}}\,\sin(4\varphi_{B})\\\
&\phantom{=}-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(4\,\varphi_{B})-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(6\,\varphi_{B})+\frac{\varepsilon^{2}\,J_{0}^{2}}{512\,\omega_{0}^{2}}\,\sin(8\,\varphi_{B})\,,\\\\[5.16663pt]
p(t)&=-\frac{3\,\varepsilon\,J_{0}^{2}}{8\,\omega_{0}}+\frac{13\,\varepsilon^{2}\,J_{0}^{3}}{16\,\omega_{0}^{2}}+\frac{\varepsilon\,J_{0}^{2}}{2\,\omega_{0}}\,\cos(2\,\varphi_{B})-\frac{33\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(2\,\varphi_{B})\\\
&\phantom{=}-\frac{\varepsilon\,J_{0}^{2}}{8\omega_{0}}\,\cos(4\,\varphi_{B})+\frac{3\,\varepsilon^{2}\,J_{0}^{3}}{16\,\omega_{0}^{2}}\,\cos(4\,\varphi_{B})+\frac{\varepsilon^{2}\,J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(6\,\varphi_{B})\,,\\\
J(t)&=J_{0}+p(t)\,.\end{split}$ (22)

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 $\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 ${\bm{m}}\cdot{\bm{\omega}_{0}}$, where
${\bm{m}}\in\mathbb{Z}^{N}$, $|{\bm{m}}|\neq 0$, and ${\bm{\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 ${\bm{\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).

(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
${{\bm{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 (20).

(iii) 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 (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

$\displaystyle{\dot{q}(t)=\frac{\partial\mathcal{H}}{\partial
J}=\omega_{0}+\frac{3\,\varepsilon\,J}{4}-\varepsilon\,J\,\cos(2q)+\frac{\varepsilon\,J\,}{4}\cos(4q)}$
(23a) $\displaystyle{\dot{J}(t)=-\frac{\partial\mathcal{H}}{\partial
q}=-\varepsilon\,J^{2}\sin(2q)+\frac{\varepsilon\,J^{2}}{2}\sin(4q)}\,.$ (23b)

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

$\omega=\omega_{0}-\varepsilon\,a_{1}(J_{0})-\varepsilon^{2}\,a_{2}(J_{0})+\ldots.$
(24)

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

$\begin{split}q(\varphi)&=q_{0}(\varphi)+\varepsilon\,q_{1}(\varphi)+\varepsilon^{2}\,q_{2}(\varphi)+\ldots\,,\\\
J(\varphi)&=J_{0}(\varphi)+\varepsilon\,J_{1}(\varphi)+\varepsilon^{2}\,J_{2}(\varphi)+\ldots\,.\end{split}$
(25)

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
${\bm{\varphi}}\in\mathbb{T}^{N}$, with ${\bm{\varphi}}={\bm{\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
${\bm{\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:

$\omega\,\frac{dq(\varphi)}{d\varphi}=\frac{dq(t)}{dt}\,,\qquad\qquad\omega\,\frac{dJ(\varphi)}{d\varphi}=\frac{dp(t)}{dt}\,.$

Substituting (24) and (25) in the equations of motion (23a) and (23b) leads to
the following expressions (up to order $2$ in $\varepsilon\,$)

$\displaystyle\left(\omega_{0}-\varepsilon\,a_{1}-\varepsilon^{2}\,a_{2}\right)\left(\frac{dq_{0}(\varphi)}{d\varphi}+\varepsilon\,\frac{dq_{1}(\varphi)}{d\varphi}+\varepsilon^{2}\,\frac{dq_{2}(\varphi)}{d\varphi}\right)$
(26a)
$\displaystyle=\omega_{0}+\frac{3\,\varepsilon\,J_{0}(\varphi)}{4}-\varepsilon\,J_{0}(\varphi)\,\cos(2\,q_{0}(\varphi))+\frac{\varepsilon\,J_{0}(\varphi)}{4}\,\cos(4\,q_{0}(\varphi))+\frac{3\,\varepsilon^{2}\,J_{1}(\varphi)}{4}-\varepsilon^{2}\,J_{1}(\varphi)\,\cos(2\,q_{0}(\varphi))$
$\displaystyle\phantom{=}+\frac{\varepsilon^{2}\,J_{1}(\varphi)}{4}\,\cos(4\,q_{0}(\varphi))+2\,\varepsilon^{2}\,J_{0}(\varphi)\,q_{1}(\varphi)\sin(2\,q_{0}(\varphi))-\varepsilon^{2}\,J_{0}(\varphi)\,q_{1}(\varphi)\,\sin(4\,\,q_{0}(\varphi))\leavevmode\nobreak\
\leavevmode\nobreak\ ,$
$\displaystyle\left(\omega_{0}-\varepsilon\,a_{1}-\varepsilon^{2}\,a_{2}\right)\left(\frac{dJ_{0}(\varphi)}{d\varphi}+\varepsilon\,\frac{dJ_{1}(\varphi)}{d\varphi}+\varepsilon^{2}\,\frac{dJ_{2}(\varphi)}{d\varphi}\right)$
(26b)
$\displaystyle=-\varepsilon\,J_{0}(\varphi)^{2}\,\sin(2\,q_{0}(\varphi))+\frac{\varepsilon\,J_{0}(\varphi)^{2}}{2}\,\sin(4\,q_{0}(\varphi))-2\,\varepsilon^{2}\,J_{0}(\varphi)\,J_{1}(\varphi)\,\sin(2\,q_{0}(\varphi))$
$\displaystyle\phantom{=}+\varepsilon^{2}\,J_{0}(\varphi)\,J_{1}(\varphi)\,\sin(4\,q_{0}(\varphi))-2\,\varepsilon^{2}\,J_{0}(\varphi)^{2}\,q_{1}(\varphi)\,\cos(2\,q_{0}(\varphi))+2\,\varepsilon^{2}\,J_{0}(\varphi)^{2}\,q_{1}(\varphi)\,\cos(4\,q_{0}(\varphi))\,.$

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

$\omega_{0}\,\frac{dq_{0}(\varphi)}{d\varphi}=\omega_{0}\,,\qquad\omega_{0}\,\frac{dJ_{0}(\varphi)}{d\varphi}=0\quad\Longrightarrow\quad
q_{0}(\varphi)=\varphi\,,\qquad J_{0}(\varphi)=J_{0}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ .$

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

$\begin{split}&\omega_{0}\,\frac{dq_{1}(\varphi)}{d\varphi}-a_{1}=\frac{3\,J_{0}}{4}-J_{0}\cos(2\varphi)+\frac{J_{0}}{4}\cos(4\varphi)\\\
&\omega_{0}\,\frac{dJ_{1}(\varphi)}{d\varphi}=-J_{0}^{2}\sin(2\varphi)+\frac{J_{0}^{2}}{2}\sin(4\,\varphi)\,.\end{split}$
(27)

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

$q_{1}(\varphi)=-\frac{J_{0}}{2\,\omega_{0}}\,\sin(2\varphi)+\frac{J_{0}}{16\,\omega_{0}}\,\sin(4\varphi)\leavevmode\nobreak\
,\qquad
J_{1}(\varphi)=-\frac{3\,J_{0}^{2}}{8\,\omega_{0}}+\frac{J_{0}^{2}}{2\,\omega_{0}}\,\cos(2\varphi)-\frac{J_{0}^{2}}{8\,\omega_{0}}\,\cos(4\varphi)\,.$

The process can be repeated at subsequent orders. We leave to the reader to
verify the result at second order, leading eventually to

$\begin{split}q(t)&=\varphi-\frac{\varepsilon\,J_{0}}{2\,\omega_{0}}\,\sin(2\varphi)+\frac{31\,\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(2\,\varphi)+\frac{\varepsilon\,J_{0}}{16\,\omega_{0}}\,\sin(4\varphi)\\\
&\phantom{=}-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(4\,\varphi)-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(6\,\varphi)+\frac{\varepsilon^{2}\,J_{0}^{2}}{512\,\omega_{0}^{2}}\,\sin(8\,\varphi)\,,\\\\[5.16663pt]
J(t)&=J_{0}-\frac{3\,\varepsilon\,J_{0}^{2}}{8\,\omega_{0}}+\frac{13\,\varepsilon^{2}\,J_{0}^{3}}{16\,\omega_{0}^{2}}+\frac{\varepsilon\,J_{0}^{2}}{2\,\omega_{0}}\,\cos(2\,\varphi)-\frac{33\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(2\,\varphi)\\\
&\phantom{=}-\frac{\varepsilon\,J_{0}^{2}}{8\omega_{0}}\,\cos(4\,\varphi)+\frac{3\,\varepsilon^{2}\,J_{0}^{3}}{16\,\omega_{0}^{2}}\,\cos(4\,\varphi)+\frac{\varepsilon^{2}\,J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(6\,\varphi)\,,\end{split}$
(28)

with $\varphi=\omega\,t$ and
$\displaystyle{\omega=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega_{0}}}\,$.
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 $\omega_{0}d/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. 222The 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$.  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}(\mathbf{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
$\mathbf{\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 (23a) and (23b) represented as a trigonometric series and
evolving according to a given pre-selected frequency $\omega$. Inverting the
expansion (24) we obtain

$\omega_{0}=\omega+\varepsilon\,a_{1}(J_{0})+\varepsilon^{2}\,a_{2}(J_{0})+\ldots$
(29)

Also, as before we expand the solution as

$\begin{split}q(t)&=q_{0}(t)+\varepsilon\,q_{1}(t)+\varepsilon^{2}\,q_{2}(t)+\ldots\\\
J(t)&=J_{0}(t)+\varepsilon\,J_{1}(t)+\varepsilon^{2}\,J_{2}(t)+\ldots\,.\end{split}$
(30)

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 $\varepsilon\,$)

$\displaystyle\dot{q_{0}}(t)+\varepsilon\,\dot{q_{1}}(t)+\varepsilon^{2}\,\dot{q_{2}}(t)=\omega+\varepsilon\,a_{1}+\frac{3\,\varepsilon\,J_{0}(t)}{4}-\varepsilon\,J_{0}(t)\,\cos(2\,q_{0}(t))+\frac{\varepsilon\,J_{0}(t)}{4}\,\cos(4q_{0}(t))$
(31a)
$\displaystyle\phantom{\dot{q_{0}}(t)+\varepsilon\,\dot{q_{1}}(t)+\varepsilon^{2}\,\dot{q_{2}}(t)=}+\frac{3\,\varepsilon^{2}\,J_{1}(t)}{4}+\varepsilon^{2}\,a_{2}-\varepsilon^{2}\,J_{1}(t)\,\cos(2\,q_{0}(t))+\frac{\varepsilon^{2}\,J_{1}(t)}{4}\,\cos(4\,q_{0}(t))$
$\displaystyle\phantom{\dot{q_{0}}(t)+\varepsilon\,\dot{q_{1}}(t)+\varepsilon^{2}\,\dot{q_{2}}(t)=}+2\,\varepsilon^{2}\,J_{0}(t)\,q_{1}(t)\,\sin(2\,q_{0}(t))-\varepsilon^{2}\,J_{0}(t)\,q_{1}(t)\,\sin(4\,q_{0}(t))\,,$
$\displaystyle\dot{J_{0}}(t)+\varepsilon\,\dot{J_{1}}(t)+\varepsilon^{2}\,\dot{J_{2}}(t)=-\varepsilon\,J_{0}(t)^{2}\,\sin(2\,q_{0}(t))+\frac{\varepsilon\,J_{0}(t)^{2}}{2}\,\sin(4q_{0}(t))$
(31b)
$\displaystyle\phantom{\dot{J_{0}}(t)+\varepsilon\,\dot{J_{1}}(t)+\varepsilon^{2}\,\dot{J_{2}}(t)=}-2\,\varepsilon^{2}\,J_{0}(t)\,J_{1}(t)\,\sin(2\,q_{0}(t))+\varepsilon^{2}\,J_{0}(t)\,J_{1}(t)\,\sin(4q_{0}(t))$
$\displaystyle\phantom{\dot{J_{0}}(t)+\varepsilon\,\dot{J_{1}}(t)+\varepsilon^{2}\,\dot{J_{2}}(t)=}-2\,\varepsilon^{2}\,J_{0}(t)^{2}\,q_{1}(t)\,\cos(2q_{0}(t))+2\,\varepsilon^{2}\,J_{0}(t)^{2}\,q_{1}(t)\,\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)=\omega\,t\,,\qquad 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 have

$\begin{split}&\displaystyle{\dot{q_{1}}(t)=a_{1}+\frac{3\,J_{0}}{4}-J_{0}\,\cos(2\,\omega\,t)+\frac{J_{0}\,}{4}\cos(4\,\omega\,t)}\,,\\\\[7.3194pt]
&\displaystyle{\dot{J_{1}}(t)=-J_{0}^{2}\sin(2\,\omega\,t)+\frac{J_{0}^{2}}{2}\sin(4\,\omega\,t)}\,\end{split}$
(32)

implying $a_{1}=-3\,J_{0}/4\,$ and

$q_{1}(t)=-\frac{J_{0}}{2\,\omega}\,\sin(2\,\omega\,t)+\frac{J_{0}}{16\,\omega}\,\sin(4\,\omega\,t)\,,\qquad
J_{1}(t)=-\frac{3\,J_{0}^{2}}{8\,\omega}+\frac{J_{0}^{2}}{2\,\omega}\,\cos(2\,\omega\,t)-\frac{J_{0}^{2}}{8\,\omega}\,\cos(4\,\omega\,t).$

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

$\begin{split}q(t)&=\omega\,t-\frac{\varepsilon\,J_{0}}{2\,\omega}\,\sin(2\,\omega\,t)+\frac{\varepsilon\,J_{0}}{16\,\omega}\,\sin(4\,\omega\,t)+\frac{19\,\varepsilon^{2}J_{0}^{2}}{32\,\omega^{2}}\,\sin(2\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{2}}{64\,\omega^{2}}\,\sin(4\,\omega\,t)\\\
&\phantom{=}-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega^{2}}\,\sin(6\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{2}}{512\,\omega^{2}}\,\sin(8\,\omega\,t)\,,\\\\[7.3194pt]
J(t)&=J_{0}-\frac{3\,\varepsilon\,J_{0}^{2}}{8\,\omega}+\frac{\varepsilon\,J_{0}^{2}}{2\,\omega}\,\cos(2\,\omega\,t)-\frac{\varepsilon\,J_{0}^{2}}{8\,\omega}\,\cos(4\,\omega\,t)+\frac{17\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\\\
&\phantom{=}-\frac{21\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(2\,\omega\,t)+\frac{3\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(4\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(6\,\omega\,t)\,,\\\\[5.16663pt]
\omega&=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega}\,.\end{split}$
(33)

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\omega t),\sin(m\omega 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, $\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
${\bm{m}}\cdot{\bm{\omega}}$. Thus, the method can proceed even when the
original frequencies of the $N$ linear oscillators ${\bm{\omega}_{0}}$ are
resonant, as long as we impose a non-resonant detuning in the adopted values
for the final frequencies ${\bm{\omega}}$ on the torus.

(iii) 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.

(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

$\omega_{0,j}-\omega_{j}=\sum_{n=1}^{\infty}\varepsilon^{n}a_{n,j}(J_{0,1},\ldots,J_{0,N})$
(34)

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 ${\bm{J}_{0}}\equiv(J_{0,1},\ldots,J_{0,N})$. Thus, the series (34) can be
formally inverted, yielding

$J_{0,j}={1\over\varepsilon}\sum_{n=1}^{\infty}P_{n,j}({\bm{\omega}}-{\bm{\omega}_{0}})$
(35)

where the functions $P_{n,j}({\bm{\omega}}-{\bm{\omega}_{0}})$ are polynomial
of degree $n$ in the differences $({\bm{\omega}}-{\bm{\omega}_{0}})$. Note
that for the inverse series to converge we must require the difference
$|({\bm{\omega}}-{\bm{\omega}_{0}})|$ to be smaller than $\varepsilon$, a fact
which limits how far we can detune ${\bm{\omega}}$ from ${\bm{\omega}_{0}}$ in
order to be able to specify the corresponding initial condition
${\bm{J}_{0}}$. An obvious choice is
$|({\bm{\omega}}-{\bm{\omega}_{0}})|=\mathcal{O}(\varepsilon^{2})$. At any
rate, we emphasize that the convergence of the inverse series (34) is an open
problem, crucial to the applications.333In the lack of a known answer to this
problem, we can adopt two different attitudes: we may fix ${\bm{J}_{0}}$ and
claim that our series answer the question of what was the value of
${\bm{\omega}_{0}}$ in the original Hamiltonian corresponding to the motion on
the particular torus with frequency ${\bm{\omega}}$ (see [20]). On the other
hand, in real world applications the only certainty we have is about our
initial model (!), i.e., ${\bm{\omega}_{0}}$, so, we must devise a method to
propose values of ${\bm{\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 ${\bm{\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 (19) with the series
(29). This yields (up to second order, apart from constants)

$\displaystyle\begin{split}\mathcal{H}(q,p)&=(\omega+\varepsilon
a_{1}+\varepsilon^{2}a_{2}+\ldots)\,p+\frac{3\,\varepsilon\,J_{0}^{2}}{8}+\frac{3\,J_{0}\,\varepsilon}{4}p+\frac{3\,\varepsilon}{8}p^{2}-\frac{\varepsilon\,J_{0}^{2}}{2}\cos(2q)-\varepsilon
J_{0}p\cos(2q)\\\
&\phantom{=}-\frac{\varepsilon}{2}p^{2}\cos(2q)+\frac{\varepsilon
J_{0}^{2}}{8}\cos(4q)+\frac{\varepsilon\,J_{0}}{4}p\cos(4q)+\frac{\varepsilon}{8}p^{2}\cos(4q)+\ldots\,.\end{split}$
(36)

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

$\displaystyle\left\langle h_{1,1}^{(0)}\right\rangle_{{q}}=\left\langle
a_{1}\,p-J_{0}\,p\,\cos(2q)+\frac{J_{0}}{4}\,p\cos(4q)+\frac{3\,J_{0}\,p}{4}\right\rangle_{{q}}=0\leavevmode\nobreak\
\leavevmode\nobreak\ ,$ (37)

leading to $a_{1}=-\frac{3\,J_{0}}{4}$. Now, we insert this expression for
$a_{1}$ in the Hamiltonian (36), leading to

$\begin{split}\mathcal{H}(q,p)&=\omega\,p+\frac{3\,\varepsilon\,J_{0}^{2}}{8}+\frac{3\,\varepsilon}{8}p^{2}-\frac{\varepsilon\,J_{0}^{2}}{2}\cos(2q)-\varepsilon
J_{0}p\cos(2q)-\frac{\varepsilon}{2}p^{2}\cos(2q)\\\
&\phantom{=}+\frac{\varepsilon
J_{0}^{2}}{8}\cos(4q)+\frac{\varepsilon\,J_{0}}{4}p\cos(4q)+\frac{\varepsilon}{8}p^{2}\cos(4q)+\varepsilon^{2}\,a_{2}\,p\,.\end{split}$
(38)

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

$h_{1,0}^{(0)}=\frac{3\,J_{0}^{2}}{8}-\frac{J_{0}^{2}}{2}\cos(2q)+\frac{J_{0}^{2}}{8}\cos(4q)\,,$

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

$X^{(1)}(q)=\frac{J_{0}^{2}}{32\,\omega}\,\sin(4q)-\frac{J_{0}^{2}}{4\,\omega}\,\sin(2q)\,.$

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

$\widehat{\mathcal{H}}^{(1)}=\exp\left(L_{\varepsilon\,\chi_{1}^{(1)}(q)}\right)\mathcal{H}\,,$

given by

$\begin{split}\widehat{\mathcal{H}}^{(1)}&=\omega\,p+\frac{3\,\varepsilon}{8}\,p^{2}-\varepsilon\,J_{0}\,p\cos(2q)-\frac{\varepsilon}{2}\,p^{2}\cos(2q)+\frac{\varepsilon\,J_{0}}{4}\,p\cos(4q)+\frac{\varepsilon}{8}\,p^{2}\cos(4q)\\\
&\phantom{=}-\frac{17\,\varepsilon^{2}\,J_{0}^{3}}{64\,\omega}+\varepsilon^{2}\,a_{2}\,p-\frac{35\,\varepsilon^{2}\,J_{0}^{2}}{64\omega}\,p+\frac{\varepsilon^{2}\,J_{0}^{3}}{2\,\omega}\,\cos(2q)+\frac{7\,\varepsilon^{2}\,J_{0}^{2}\,p}{8\,\omega}\,\cos(2q)-\frac{11\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega}\,\cos(4q)\\\
&\phantom{=}-\frac{7\,\varepsilon^{2}\,J_{0}^{2}\,p}{16\,\omega}\,\cos(4q)+\frac{\varepsilon^{2}\,J_{0}^{3}}{8\,\omega}\,\cos(6q)+\frac{\varepsilon^{2}\,J_{0}^{2}\,p}{8\,\omega}\,\cos(6q)-\frac{\varepsilon^{2}\,J_{0}^{3}}{64\,\omega}\,\cos(8q)-\frac{\varepsilon^{2}\,J_{0}^{2}\,p}{64\,\omega}\,\cos(8q)\,.\end{split}$
(39)

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

$\widehat{h}_{1,1}^{(1)}=-J_{0}\,p\cos(2q)+\frac{1}{4}\,J_{0}\,p\cos(4q)\,$

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

$\tilde{\chi}_{2}^{(1)}(q,p)=-\frac{J_{0}\,p}{2\,\omega}\,\sin(2q)+\frac{J_{0}\,p}{16\,\omega}\,\sin(4q)\,.$

The new Hamiltonian is

$\mathcal{H}^{(1)}=\exp\left(L_{\varepsilon\,\chi_{2}^{(1)}(q,p)}\right)\widehat{\mathcal{H}}^{(1)}\,$

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

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

$\displaystyle a_{2}$ $\displaystyle=\frac{69\,J_{0}^{2}}{64\,\omega}\,,$ (40)
$\displaystyle\chi_{1}^{(2)}$
$\displaystyle=-\frac{17\,J_{0}^{3}\,q}{64\,\omega^{2}}+\frac{J_{0}^{3}}{4\,\omega^{2}}\,\sin(2q)-\frac{11\,J_{0}^{3}}{128\,\omega^{2}}\,\sin(4q)+\frac{J_{0}^{3}}{48\,\omega^{2}}\,\sin(6q)-\frac{J_{0}^{3}}{512\,\omega^{2}}\,\sin(8q)\,,$
$\displaystyle\chi_{2}^{(2)}$
$\displaystyle=\frac{37\,J_{0}^{2}\,p}{64\,\omega^{2}}\,\sin(2q)-\frac{7\,J_{0}^{2}\,p}{64\,\omega^{2}}\,\sin(4q)+\frac{J_{0}^{2}\,p}{64\,\omega^{2}}\,\sin(6q)-\frac{J_{0}^{2}\,p}{512\,\omega^{2}}\,\sin(8q)\,,$

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

$\displaystyle Z^{(2)}(\tilde{q},\tilde{p})$
$\displaystyle=\omega\,\tilde{p}+\frac{3\,\varepsilon\,\tilde{p}^{2}}{8}-\frac{\varepsilon\,\tilde{p}^{2}}{2}\,\cos(2\tilde{q})+\frac{\varepsilon\,\tilde{p}^{2}}{8}\,\cos(4\tilde{q})-\frac{51\,\varepsilon^{2}\,J_{0}\,\tilde{p}^{2}}{64\,\omega}+\frac{37\,\varepsilon^{2}\,J_{0}\,\tilde{p}^{2}}{32\,\omega}\,\cos(2\tilde{q})$
$\displaystyle\phantom{=}-\frac{7\,\varepsilon^{2}\,J_{0}\,\tilde{p}^{2}}{16\,\omega}\,\cos(4\tilde{q})+\frac{3\,\varepsilon^{2}\,J_{0}\,\tilde{p}^{2}}{32\,\omega}\,\cos(6\tilde{q})-\frac{\varepsilon^{2}\,J_{0}\,\tilde{p}^{2}}{64\,\omega}\,\cos(8\tilde{q})\,$

($(\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

$Z^{(2)}(\tilde{q},\tilde{p})=\omega\,\tilde{p}+\varepsilon\,R_{1}(\tilde{q},\tilde{p})+\varepsilon^{2}\,R_{2}(\tilde{q},\tilde{p})\,,$

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

$\begin{cases}&\displaystyle{\dot{\tilde{q}}=\frac{\partial
Z^{(2)}}{\partial\tilde{p}}=\omega+\varepsilon\,\frac{\partial
R_{1}(\tilde{q},\tilde{p})}{\partial\tilde{p}}+\varepsilon^{2}\,\frac{\partial
R_{2}(\tilde{q},\tilde{p})}{\partial\tilde{p}}=\omega+\mathcal{O}(||\,\tilde{p}\,||)}\\\
&\displaystyle{\dot{\tilde{p}}=-\frac{\partial
Z^{(2)}}{\partial\tilde{q}}=-\varepsilon\,\frac{\partial
R_{1}(\tilde{q},\tilde{p})}{\partial\tilde{q}}-\varepsilon^{2}\,\frac{\partial
R_{2}(\tilde{q},\tilde{p})}{\partial\tilde{q}}=\mathcal{O}(||\,\tilde{p}\,||^{2})}\end{cases}\,.$

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

$\displaystyle
q=\mathrm{exp}L_{\varepsilon^{2}\,\chi_{2}^{(2)}}\mathrm{exp}L_{\varepsilon\,\chi_{2}^{(1)}}\,\tilde{q}\,,$
$\displaystyle
p=\mathrm{exp}L_{\varepsilon^{2}\,\chi_{2}^{(2)}}\mathrm{exp}L_{\varepsilon^{2}\,\chi_{1}^{(2)}}\mathrm{exp}L_{\varepsilon\,\chi_{2}^{(1)}}\mathrm{exp}L_{\varepsilon\,\chi_{1}^{(1)}}\,\tilde{p}\,.$
(41)

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:

$\begin{split}q(t)&=\omega\,t-\frac{\varepsilon\,J_{0}}{2\,\omega}\,\sin(2\,\omega\,t)+\frac{\varepsilon\,J_{0}}{16\,\omega}\,\sin(4\,\omega\,t)+\frac{19\,\varepsilon^{2}J_{0}^{2}}{32\,\omega^{2}}\,\sin(2\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{2}}{64\,\omega^{2}}\,\sin(4\,\omega\,t)\\\
&\phantom{=}-\frac{\varepsilon^{2}\,J_{0}^{2}}{32\,\omega^{2}}\,\sin(6\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{2}}{512\,\omega^{2}}\,\sin(8\,\omega\,t)\,,\\\\[7.3194pt]
p(t)&=-\frac{3\,\varepsilon\,J_{0}^{2}}{8\,\omega}+\frac{\varepsilon\,J_{0}^{2}}{2\,\omega}\,\cos(2\,\omega\,t)-\frac{\varepsilon\,J_{0}^{2}}{8\,\omega}\,\cos(4\,\omega\,t)+\frac{17\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\\\
&\phantom{=}-\frac{21\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(2\,\omega\,t)+\frac{3\,\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(4\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(6\,\omega\,t)\,,\\\
J(t)&=J_{0}+p(t)\,.\end{split}$ (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

$\omega=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega}\,.$

These expressions are identical to the ones found by the Lindstedt method of
Subsection 2.3 (see (33)).

###  2.5 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
$\varepsilon^{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
$\varepsilon^{2}\,$.

Method: | $\mathcal{O}(1)$ | $\mathcal{O}(\varepsilon)$ | $\mathcal{O}(\varepsilon^{2})$
---|---|---|---
(i) ‘Analogous to Birkhoff’ | $\displaystyle{\omega\,t}$ | $\displaystyle{-\frac{J_{0}}{2\,\omega_{0}}\,\sin(2\,\omega\,t)+\frac{J_{0}}{16\,\omega_{0}}\,\sin(4\,\omega\,t)}$ | $\displaystyle{\frac{31\,J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(2\,\omega\,t)-\frac{J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(4\,\omega\,t)\quad}\atop{\displaystyle{-\frac{J_{0}^{2}}{32\,\omega_{0}^{2}}\,\sin(6\,\omega\,t)+\frac{J_{0}^{2}}{512\,\omega_{0}^{2}}\,\sin(8\,\omega\,t)}}$
(ii) ‘Analogous to Kolmogorov’ | $\displaystyle{\omega\,t}$ | $\displaystyle{-\frac{J_{0}}{2\,\omega}\,\sin(2\,\omega\,t)+\frac{J_{0}}{16\,\omega}\,\sin(4\,\omega\,t)}$ | $\displaystyle{\frac{19\,J_{0}^{2}}{32\,\omega^{2}}\,\sin(2\,\omega\,t)+\frac{J_{0}^{2}}{64\,\omega^{2}}\,\sin(4\,\omega\,t)\quad}\atop{\displaystyle{-\frac{J_{0}^{2}}{32\,\omega^{2}}\,\sin(6\,\omega\,t)+\frac{J_{0}^{2}}{512\,\omega^{2}}\,\sin(8\,\omega\,t)}}$

Table 2: Comparison between the series terms for the solution $J(t)\,$ in the
Lindstedt series obtained by the methods (i) and (ii) up to order
$\varepsilon^{2}\,$.

Method: | $\mathcal{O}(1)$ | $\mathcal{O}(\varepsilon)$ | $\mathcal{O}(\varepsilon^{2})$
---|---|---|---
(i) ‘Analogous to Birkhoff’ | $\displaystyle{J_{0}}$ | $\displaystyle{-\frac{3\,J_{0}^{2}}{8\,\omega_{0}}+\frac{J_{0}^{2}}{2\,\omega_{0}}\,\cos(2\,\omega\,t)-\frac{J_{0}^{2}}{8\omega_{0}}\,\cos(4\,\omega\,t)}$ | $\displaystyle{\frac{13\,J_{0}^{3}}{16\,\omega_{0}^{2}}-\frac{33\,J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(2\,\omega\,t)\qquad\qquad\quad}\atop{\displaystyle{+\frac{3\,J_{0}^{3}}{16\,\omega_{0}^{2}}\,\cos(4\,\omega\,t)+\frac{J_{0}^{3}}{32\,\omega_{0}^{2}}\,\cos(6\,\omega\,t)}}$
(ii)‘Analogous to Kolmogorov’ | $\displaystyle{J_{0}}$ | $\displaystyle{-\frac{3\,J_{0}^{2}}{8\,\omega}+\frac{J_{0}^{2}}{2\,\omega}\,\cos(2\,\omega\,t)-\frac{J_{0}^{2}}{8\,\omega}\,\cos(4\,\omega\,t)}$ | $\displaystyle{\frac{17\,J_{0}^{3}}{32\,\omega^{2}}-\frac{21\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(2\,\omega\,t)\qquad\qquad\quad}\atop{\displaystyle{+\frac{3\,J_{0}^{3}}{32\,\omega^{2}}\,\cos(4\,\omega\,t)+\frac{J_{0}^{3}}{32\,\omega^{2}}\,\cos(6\,\omega\,t)}}$

From 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
$\mathcal{O}(\varepsilon)\,$, except for the fact that in method (i) 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:

$\displaystyle\omega=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega_{0}}\,,$
$\displaystyle\omega=\omega_{0}+\frac{3\,\varepsilon\,J_{0}}{4}-\frac{69\,\varepsilon^{2}\,J_{0}^{2}}{64\,\omega}\,,$

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\geq 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 $\omega_{0}$ and $\omega\,$, we reverse the series (24) according to

$d\omega=-\sum_{i\geq
0}\varepsilon^{i}\,a_{i}(J_{0}):=f(J_{0})\quad\Longrightarrow\quad
J_{0}=f^{-1}(d\omega)\,,$ (43)

where $d\omega=\omega-\omega_{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 $\varepsilon^{4}\,$. We analyze three different cases:

  1. 1)

$\varepsilon=1\,$, $\omega_{0}=1$ and $\omega=1.002\,$;

  2. 2)

$\varepsilon=1\,$, $\omega_{0}=1$ and $\omega=1.02\,$;

  3. 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 (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
$\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
(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 $\mathrm{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\omega$ is of order $10^{-3}$ and raising up to two orders of magnitude when
$d\omega$ becomes of the order of unity.

(a)

(b)

(c)

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)
$\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\,$.

###  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

$\mathcal{H}({\bm{q}},{\bm{J}})={\bm{\omega}_{0}}\cdot{\bm{J}}+\varepsilon
h({\bm{q}},{\bm{J}};\varepsilon)$

where the development in powers of the variables ${\bm{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

$\mathcal{H}({\bm{q}},{\bm{J}})\,=\,{\bm{\omega}}_{0}\cdot{\bm{J}}+\sum_{i\geq
1}\sum_{s\geq\,3}\,\varepsilon^{i}\tilde{h}_{i,\,s}({\bm{q}},{\bm{J}})$ (44)

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

$\begin{split}\mathcal{H}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,&=\,{\bm{\omega}}_{0}\cdot{\bm{J}_{0}}+{\bm{\omega}}_{0}\cdot{\bm{p}}+\sum_{i\geq
1}\sum_{s\geq\,3}\,\varepsilon^{i}\tilde{h}_{i,\,s}({\bm{q}},{\bm{J}_{0}}+{\bm{p}})\\\
&=\,{\bm{\omega}}_{0}\cdot{\bm{J}_{0}}+{\bm{\omega}}_{0}\cdot{\bm{p}}+\sum_{i\geq
1}\sum_{s\geq\,3}\sum_{k=0}^{n}\,\varepsilon^{i}\frac{\tilde{h}^{(k)}_{i,\,s}({\bm{q}},{\bm{J}_{0}}+{\bm{p}})|_{{\bm{p}}={\bm{0}}}}{k!}\,{\bm{p}}^{k}\,,\end{split}$
(45)

where $\tilde{h}^{(k)}$ are the $k-$th derivatives of $\tilde{h}$ with respect
to ${\bm{p}}={\bm{J}}-{\bm{J}_{0}}\,$.

Apart from constants, the Hamiltonian can be written in the compact notation

$\mathcal{H}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,=\,{\bm{\omega}}_{0}\cdot{\bm{p}}+\sum_{i\geq\,1}\,\varepsilon^{i}\,h_{i}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,.$
(46)

The algorithm allows to compute quasi-periodic orbits with a frequency
${\bm{\omega}}$ fixed in advance, given by

${\bm{\omega}}={\bm{\omega}}_{0}-\sum_{i\geq\,1}\,\varepsilon^{i}\,{\bm{a}}_{i}({\bm{J}_{0}})\,,$
(47)

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

Substituting the series (47) in the Hamiltonian (46) we arrive at:

$\mathcal{H}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,=\,{\bm{\omega}}\cdot{\bm{p}}+\sum_{i\geq\,1}\,\varepsilon^{i}\,h_{i}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\sum_{i\geq\,1}\,\varepsilon^{i}\,{\bm{a}}_{i}({\bm{J}_{0}})\cdot{\bm{p}}\,.$
(48)

In the following we use the notation

$\varepsilon^{i}\,h_{i}^{(k)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})=\varepsilon^{i}\left(h_{i,0}^{(k)}({\bm{q}};{\bm{J}_{0}})+h_{i,1}^{(k)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\sum_{j\geq\,2}h_{i,j}^{(k)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\right)\,,$
(49)

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

$\begin{split}&\phantom{h_{i,j}^{(k)}\,:}\,\,\>i=\textit{degree of the
corresponding $\varepsilon$}\,,\\\ &h_{i,j}^{(k)}\,:\>j=\textit{degree of
${\bm{p}}$}\,,\\\ &\phantom{h_{i,j}^{(k)}\,:}\,\>k=\textit{step of the
algorithm}\,.\end{split}$ (50)

We now have the following

######  Proposition 1.

Assume the vector ${\bm{\omega}}$ is non-resonant. There exist Lie generating
functions $\chi_{1}^{(r)}$, $\chi_{2}^{(r)}$ such that, after $r$
normalization steps, the Hamiltonian (48) is given by the formal series:

$\displaystyle\mathcal{H}^{(r)}({\bm{q}}^{(r)},{\bm{p}}^{(r)};{\bm{J}}_{0})$
$\displaystyle=\bigg{(}\exp(L_{\varepsilon^{r}\chi_{2}^{(r)}})\big{(}\exp(L_{\varepsilon^{r}\chi_{1}^{(r)}})\mathcal{H}^{(r-1)}\big{)}\bigg{)}\bigg{|}_{\begin{subarray}{c}{\bm{q}^{(r-1)}}={\bm{q}}^{(r)}\\\
{\bm{p}^{(r-1)}}={\bm{p}}^{(r)}\end{subarray}}\,$
$\displaystyle=\,{\bm{\omega}}\cdot{\bm{p}}+\varepsilon\,C_{1}+\ldots+\varepsilon^{r-1}\,C_{r-1}+\varepsilon^{r}\,C_{r}$
$\displaystyle\phantom{=}+\varepsilon\,h_{1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots+\varepsilon^{r}\,h_{r}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})$
$\displaystyle\phantom{=}+\varepsilon^{r+1}\,h_{r+1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\varepsilon^{r+1}\,{\bm{a}}_{r+1}({\bm{J}_{0}})\cdot{\bm{p}}+\ldots\,,$

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}(||{\bm{p}}||^{2})\,\forall\,j=1,\ldots,r$,
with

$\displaystyle h_{k,0}^{(r)}=C_{k}=\left\langle
h_{k,0}^{(k-1)}\right\rangle-{\bm{\omega}}\cdot{\bm{K}}^{(k)}\,$
$\displaystyle\;\;\forall\;1\leq k\leq r\,,$ $\displaystyle
h_{k,0}^{(r)}=h_{k,0}^{(r-1)}$ $\displaystyle\;\;\forall\;r+1\leq k\leq
2r-1\,,$ $\displaystyle
h_{2r,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{2r-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{1}^{(r)}({\bm{q}})}^{j}\left(h_{2r-jr,j}^{(r-1)}\right)$
$\displaystyle
h_{k,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-2}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,s}^{(r-1)}\right)$
$\displaystyle\qquad+\frac{1}{\left(\lfloor\frac{k-1}{r}\rfloor-1\right)!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor-1}\left(h_{k-\left(\lfloor\frac{k-1}{r}\rfloor-1\right)r,0}^{(r-1)}\right)$
$\displaystyle\;\;\forall\;k\geq 2r+1\,,k\neq mr\,,m\in\mathbb{N}\,,$
$\displaystyle
h_{k,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,s}^{(r-1)}\right)$ $\displaystyle\;\;\forall\;k\geq
2r+1\,,k=mr\,,m\in\mathbb{N}\,,$ $\displaystyle h_{k,1}^{(r)}=0$
$\displaystyle\;\;\forall\;1\leq k\leq r\,,$ $\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,1+s}^{(r-1)}\right)$ $\displaystyle\;\;\forall\;k\geq r+1\,,k\neq
m\,r\,,\;m\in\mathbb{N}\,,$ $\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-2}{r}\rfloor}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,1+s}^{(r-1)}\right)$
$\displaystyle\quad\phantom{h_{k,i}^{(r)}=}+\frac{m-1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(h_{r,1}^{(r-1)}\right)$
$\displaystyle\;\;\forall\;k\geq r+1\,,k=m\,r\,,\;m\in\mathbb{N}\,,$
$\displaystyle h_{k,i}^{(r)}=h_{k,i}^{(r-1)}$ $\displaystyle\;\;\forall\;1\leq
k\leq r\,,\;i\geq 2\,,$ $\displaystyle
h_{k,i}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,i+s}^{(r-1)}\right)$ $\displaystyle\;\;\forall\;k\geq r+1\,,\;i\geq 2\,.$

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

$\displaystyle\chi_{1}^{(r)}({\bm{q}})=X^{(r)}({\bm{q}})+{\bm{K}}^{(r)}\cdot{\bm{q}}\,=\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left(\frac{c^{(r-1)}_{{\bm{k}}}}{i{\bm{k}}\cdot{\bm{\omega}}}\,e^{i{\bm{k}}\cdot{\bm{q}}}-\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,{\bm{k}}\cdot{\bm{q}}\right)\,,$
$\displaystyle\chi_{2}^{(r)}({\bm{q}},{\bm{p}})=\tilde{\chi}_{2}^{(r)}({\bm{q}},{\bm{p}})+{\bm{S}}^{(r)}\cdot{\bm{p}}$
$\displaystyle\phantom{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}={\bm{Y}}^{(r)}({\bm{q}})\cdot{\bm{p}}+{\bm{S}}^{(r)}\cdot{\bm{p}}=\sum_{|{\bm{j}}|=1}\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left(\frac{c_{{\bm{j}}{\bm{k}}}^{(r)}}{i{\bm{k}}\cdot{\bm{\omega}}}\,e^{i{\bm{k}}\cdot{\bm{q}}}+\frac{i\,c_{{\bm{j}}{\bm{k}}}^{(r)}}{{\bm{k}}\cdot{\bm{\omega}}}\right){\bm{p}}^{{\bm{j}}}\,.$

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:

$\begin{split}\mathcal{H}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,=\,&{\bm{\omega}}\cdot{\bm{p}}+\varepsilon\,C_{1}+\ldots+\varepsilon^{r-1}\,C_{r-1}\\\
&+\varepsilon\,h_{1}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots+\varepsilon^{r-1}\,h_{r-1}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\\\
&+\varepsilon^{r}\,h_{r}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\varepsilon^{r}\,{\bm{a}}_{r}({\bm{J}_{0}})\cdot{\bm{p}}\\\
&+\varepsilon^{r+1}\,h_{r+1}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\varepsilon^{r+1}\,{\bm{a}}_{r+1}({\bm{J}_{0}})\cdot{\bm{p}}+\ldots\,,\end{split}$
(51)

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

$h_{i,1}^{(k)}\leftarrow
h_{i,1}^{(k)}+{\bm{a}}_{i}\cdot{\bm{p}}=\left(\nabla_{{\bm{p}}}h_{i,1}^{(k)}+{\bm{a}}_{i}\right)\cdot{\bm{p}}\quad\forall\,i\geq
r\,.$ (52)

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

$\widehat{\mathcal{H}}^{(r)}=\exp(L_{\varepsilon^{r}\,\chi_{1}^{(r)}({\bm{q}})})\mathcal{H}^{(r-1)}=\sum_{j\geq
0}\frac{(\varepsilon^{r})^{j}}{j!}L_{\chi_{1}^{(r)}({\bm{q}})}^{j}\mathcal{H}^{(r-1)}\,,$
(53)

where the generating function $\chi_{1}^{(r)}$ is defined as
$\chi_{1}^{(r)}({\bm{q}})=X^{(r)}({\bm{q}})+{\bm{K}}^{(r)}\cdot{\bm{q}}\,$,
with ${\bm{K}}^{(r)}$ an appropriate constant vector defined below, and
$X^{(r)}$ is defined through the homological equation:

$\big{\\{}{\bm{\omega}}\cdot{\bm{p}}\,,\,X^{(r)}({\bm{q}})\big{\\}}+h_{r,0}^{(r-1)}({\bm{q}})=\left\langle
h_{r,0}^{(r-1)}({\bm{q}})\right\rangle\,,$ (54)

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

$X^{(r)}({\bm{q}})=\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{c^{(r-1)}_{{\bm{k}}}}{i{\bm{k}}\cdot{\bm{\omega}}}\,e^{i{\bm{k}}\cdot{\bm{q}}}.\,$
(55)

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

$\left\langle
h_{r,1}^{(r-1)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\right\rangle=0\,.$ (56)

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

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

$X^{(r)}({\bm{q}})=\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[-\frac{i\,c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\cos({\bm{k}}\cdot{\bm{q}})+\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\sin({\bm{k}}\cdot{\bm{q}})\right]\,,$
(57)

the poisson bracket of $\chi_{1}^{(r)}({\bm{q}})$ with the i-th component
$p_{i}\,$ ($1\leq\,i\,\leq n\,$) of the vector ${\bm{p}}$ yields the following
expressions for the i-th components $K_{i}^{(r)}$ and $k_{i}$ of the vectors
${\bm{K}}^{(r)}$ and ${\bm{k}}\,$:

$\begin{split}\\{p_{i}\,,X^{(r)}({\bm{q}})\\}+\\{p_{i}\,,{\bm{K}}^{(r)}\cdot{\bm{q}}\\}&=\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[-\big{\\{}p_{i}\,,\,\frac{i\,c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\cos({\bm{k}}\cdot{\bm{q}})\big{\\}}+\big{\\{}p_{i}\,,\,\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\sin({\bm{k}}\cdot{\bm{q}})\big{\\}}\right]-K^{(r)}_{i}\\\
&=-\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[\,\frac{i\,c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,k_{i}\,\sin({\bm{k}}\cdot{\bm{q}})+\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,k_{i}\,\cos({\bm{k}}\cdot{\bm{q}})\,\right]-K^{(r)}_{i}=\\\
&=-\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{i\,c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,k_{i}\,\sin({\bm{k}}\cdot{\bm{q}})-\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,k_{i}\left(\cos({\bm{k}}\cdot{\bm{q}})-1\right)\,.\end{split}$
(58)

Hence

$K^{(r)}_{i}=-\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,k_{i}\,,$
(59)

i.e.,

${\bm{K}}^{(r)}=-\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{c^{(r-1)}_{{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,{\bm{k}}\,.$
(60)

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

$\begin{split}\widehat{\mathcal{H}}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,=\,&{\bm{\omega}}\cdot{\bm{p}}+\varepsilon\,C_{1}+\ldots+\varepsilon^{r-1}\,C_{r-1}+\varepsilon^{r}\,C_{r}\\\
&+\varepsilon\,\widehat{h}_{1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots+\varepsilon^{r-1}\,\widehat{h}_{r-1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\\\
&+\varepsilon^{r}\,\widehat{h}_{r}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\varepsilon^{r+1}\,\widehat{h}_{r+1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots\,,\end{split}$
(61)

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}(||\,{\bm{p}}\,||^{2})\,$
$\forall\,j=1,\ldots,r-1\,$, where:

$\displaystyle\widehat{h}_{k,0}^{(r)}=C_{k}=\left\langle
h_{k,0}^{(k-1)}\right\rangle-{\bm{\omega}}\cdot{\bm{K}}^{(k)}$
$\displaystyle\forall\,1\leq k\leq r\,,$ (62)
$\displaystyle\widehat{h}_{k,0}^{(r)}=h_{k,0}^{(r-1)}$
$\displaystyle\forall\>r+1\leq k\leq 2r-1\,,$ (63)
$\displaystyle\widehat{h}_{k,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{1}^{(r)}({\bm{q}})}^{j}\left(h_{k-jr,j}^{(r-1)}\right)$
$\displaystyle\forall\,k\geq 2r\,$ (64)
$\displaystyle\widehat{h}_{k,1}^{(r)}=0$ $\displaystyle\forall\>1\leq k\leq
r-1\,,$ (65) $\displaystyle\widehat{h}_{r,1}^{(r)}=h_{r,1}^{(r-1)}\,,$ (66)
$\displaystyle\widehat{h}_{k,i}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{1}^{(r)}({\bm{q}})}^{j}\left(h_{k-jr,i+j}^{(r-1)}\right)$
$\displaystyle\forall\>k\geq r+1\,,\>i\geq 1\,,$ (67)
$\displaystyle\widehat{h}_{k,i}^{(r)}=h_{k,i}^{(r-1)}$
$\displaystyle\forall\>1\leq k\leq r\,,\>i\geq 2\,.$ (68)

  * •

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

$\mathcal{H}^{(r)}=\exp(L_{\varepsilon^{r}\,\chi_{2}^{(r)}({\bm{q}},{\bm{p}})})\widehat{\mathcal{H}}^{(r)}=\sum_{j\geq
0}\frac{(\varepsilon^{r})^{j}}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\widehat{\mathcal{H}}^{(r)}\,,$
(69)

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

$\big{\\{}{\bm{\omega}}\cdot{\bm{p}}\,,\,\tilde{\chi}_{2}^{(r)}({\bm{q}},{\bm{p}})\big{\\}}+\widehat{h}_{r,1}^{(r)}({\bm{q}},{\bm{p}})=0\,.$
(70)

Setting
$\displaystyle{\widehat{h}_{r,1}^{(r)}({\bm{q}},{\bm{p}})=\sum_{|{\bm{j}}|=1}\sum_{{\bm{k}}\in\mathbb{Z}^{n}}c_{{\bm{j}}{\bm{k}}}^{(r)}{\bm{p}}^{{\bm{j}}}e^{i{\bm{k}}\cdot{\bm{q}}}}$
the solution of (70) is

$\displaystyle{\tilde{\chi}_{2}^{(r)}({\bm{q}},{\bm{p}})=\sum_{|{\bm{j}}|=1}\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{c_{{\bm{j}}{\bm{k}}}^{(r)}}{i{\bm{k}}\cdot{\bm{\omega}}}{\bm{p}}^{{\bm{j}}}\,e^{i{\bm{k}}\cdot{\bm{q}}}}.\,$
(71)

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

$\begin{split}\\{q_{i}\,,\chi_{2}^{(r)}\\}&=\\{q_{i}\,,\tilde{\chi}_{2}^{(r)}({\bm{q}},{\bm{p}})\\}+\\{q_{i}\,,{\bm{S}}^{(r)}\cdot{\bm{p}}\\}\\\
&=\sum_{|{\bm{j}}|=1}\left(\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[\big{\\{}q_{i}\,,\,\frac{c_{{\bm{j}}{\bm{k}}}^{(r)}}{i\,{\bm{k}}\cdot{\bm{\omega}}}\,{\bm{p}}^{{\bm{j}}}\,\cos({\bm{k}}\cdot{\bm{q}})\big{\\}}+\big{\\{}q_{i}\,,\,\frac{c^{(r)}_{{\bm{j}}{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,{\bm{p}}^{{\bm{j}}}\sin({\bm{k}}\cdot{\bm{q}})\big{\\}}\right]+\big{\\{}q_{i}\,\,,\,S_{{\bm{j}}}^{(r)}{\bm{p}}^{{\bm{j}}}\big{\\}}\right)\\\
&=\sum_{|{\bm{j}}|=1}\left(\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[\frac{c_{{\bm{j}}{\bm{k}}}^{(r)}}{i\,{\bm{k}}\cdot{\bm{\omega}}}\,j_{i}\,p_{i}^{j_{i}-1}\,\cos({\bm{k}}\cdot{\bm{q}})+\frac{c^{(r)}_{{\bm{j}}{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,j_{i}\,p_{i}^{j_{i}-1}\sin({\bm{k}}\cdot{\bm{q}})\right]+{S_{{\bm{j}}}^{(r)}\,j_{i}\,p_{i}^{j_{i}-1}}\right)\\\
&=\sum_{|{\bm{j}}|=1}\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\left[\frac{c_{{\bm{j}}{\bm{k}}}^{(r)}}{i\,{\bm{k}}\cdot{\bm{\omega}}}\,j_{i}\,p_{i}^{j_{i}-1}\,\left(\cos({\bm{k}}\cdot{\bm{q}})-1\right)+\frac{c^{(r)}_{{\bm{j}}{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,j_{i}\,p_{i}^{j_{i}-1}\sin({\bm{k}}\cdot{\bm{q}})\right]\,\end{split}$
(72)

we obtain

$S^{(r)}_{\bm{j}}=\sum_{{\bm{k}}\in\mathbb{Z}^{n}\setminus\\{{\bm{0}}\\}}\frac{i\,c^{(r)}_{{\bm{j}}{\bm{k}}}}{{\bm{k}}\cdot{\bm{\omega}}}\,.$
(73)

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

$\begin{split}\mathcal{H}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\,=\,&{\bm{\omega}}\cdot{\bm{p}}+\varepsilon\,C_{1}+\ldots+\varepsilon^{r-1}\,C_{r-1}+\varepsilon^{r}\,C_{r}\\\
&+\varepsilon\,h_{1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots+\varepsilon^{r}\,h_{r}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})\\\
&+\varepsilon^{r+1}\,h_{r+1}^{(r)}({\bm{q}},{\bm{p}};{\bm{J}_{0}})+\ldots\,,\end{split}$
(74)

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}(||\,{\bm{p}}\,||^{2})\,\forall\,j=1,\ldots,r\,$,
where:

$\displaystyle h_{k,0}^{(r)}=\widehat{h}_{k,0}^{(r)}$
$\displaystyle\;\;\forall\;1\leq k\leq 2r\,,$ (75) $\displaystyle
h_{k,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,0}^{(r)}\right)$
$\displaystyle\;\;\forall\;k\geq 2r+1\,,$ (76) $\displaystyle h_{k,1}^{(r)}=0$
$\displaystyle\;\;\forall\;1\leq k\leq r\,,$ (77) $\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,1}^{(r)}\right)$
$\displaystyle\;\;\forall\;k\geq r+1\,,k\neq m\,r\,(m\in\mathbb{N})\,,$ (78)
$\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,1}^{(r)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)$
$\displaystyle\;\;\forall\;k\geq r+1\,,k=m\,r\,(m\in\mathbb{N})\,.$ (79)
$\displaystyle h_{k,i}^{(r)}=\widehat{h}_{k,i}^{(r)}$
$\displaystyle\;\;\forall\;1\leq k\leq r\,,\>i\geq 2\,,$ (80) $\displaystyle
h_{k,i}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,i}^{(r)}\right)$
$\displaystyle\;\;\forall\;k\geq r+1\,,\>i\geq 2\,.$ (81)

Equation (79), using homological equation (70), can be written equivalently
as:

$\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{m-2}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,1}^{(r)}\right)+\frac{m-1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(\widehat{h}_{r,1}^{(r-1)}\right)$
$\displaystyle\;\;\forall\;k\geq r+1\,,k=m\,r\,(m\in\mathbb{N})\,$ (82)

where, from (79) we have

$\displaystyle\sum_{j=0}^{m-1}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{mr-
jr,1}^{(r)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)=$
$\displaystyle\sum_{j=0}^{m-2}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{mr-
jr,1}^{(r)}\right)+\frac{1}{(m-1)!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(\widehat{h}_{r,1}^{(r)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)=$
$\displaystyle\sum_{j=0}^{m-2}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{mr-
jr,1}^{(r)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(\widehat{h}_{r,1}^{(r)}+L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}\left({\bm{\omega}}\cdot{\bm{p}}\right)\right)+\frac{m-1}{m!}\,L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(\widehat{h}_{r,1}^{(r)}\right).$

_Second part of the proof :_ using the formulas (62)–(68) and (75)–(81), 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 (75) we have

$\displaystyle h_{k,0}^{(r)}=\widehat{h}_{k,0}^{(r)}=C_{k}=\left\langle
h_{k,0}^{(k-1)}\right\rangle-{\bm{\omega}}\cdot{\bm{K}}^{(k)}\,\qquad\leavevmode\nobreak\
\forall\;1\leq k\leq r\,,$ (83) $\displaystyle
h_{k,0}^{(r)}=\widehat{h}_{k,0}^{(r)}=h_{k,0}^{(r-1)}\qquad\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \forall\>r+1\leq k\leq 2r-1\,,$ (84)
$\displaystyle
h_{2r,0}^{(r)}=\widehat{h}_{2r,0}^{(r)}=\sum_{j=0}^{\lfloor\frac{2r-1}{r}\rfloor}\frac{1}{j!}L_{\chi_{1}^{(r)}({\bm{q}})}^{j}\left(h_{2r-jr,j}^{(r-1)}\right)=h_{2r,0}^{(r-1)}+L_{\chi_{1}^{(r)}({\bm{q}})}\left(h_{r,1}^{(r-1)}\right)\,,$
(85)

where we used Eqs (62), (63) and (64), respectively. Using (68) we have

$h_{k,i}^{(r)}=\widehat{h}_{k,i}^{(r)}=h_{k,i}^{(r-1)}\quad\forall\>1\leq
k\leq r\,,\>i\geq 2\,.$ (86)

Equation (81) can now be written in the form

$h_{k,i}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,i}^{(r)}\right)+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(\widehat{h}_{k-\lfloor\frac{k-1}{r}\rfloor
r,i}^{(r)}\right)$ (87)

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

$k-jr\geq k-\lfloor\frac{k-1}{r}\rfloor r+r\geq r+1\,,$ (88)

where we have used the inequalities $j\leq\lfloor(k-1)/r\rfloor-1$ and
$\lfloor k-1\rfloor\leq 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\leq f\leq r-1\,$ and $m\in\mathbb{N}\,$. Thus

$\displaystyle k-\lfloor\frac{k-1}{r}\rfloor
r=mr+f-\lfloor\frac{mr+f-1}{r}\rfloor r=f-\lfloor\frac{f-1}{r}\rfloor r\,.$
(89)

Then

$\displaystyle k-\lfloor\frac{k-1}{r}\rfloor r=f-\lfloor\frac{f-1}{r}\rfloor
r=r\quad\mathrm{if}\;\,f=0\,$ (90)

and

$\displaystyle 1\leq k-\lfloor\frac{k-1}{r}\rfloor
r=f-\lfloor\frac{f-1}{r}\rfloor r\leq r-1\quad\mathrm{if}\;\,1\leq f\leq
r-1\,.$ (91)

We can conclude that $1\leq k-\lfloor(k-1)/r\rfloor r\leq r\,$, implying that
the definition (68) holds for the second term in (87). We can then write (87)
as444Observe that, due to (90) and (91), 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.

$\displaystyle h_{k,i}^{(r)}$
$\displaystyle=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,i+s}^{(r-1)}\right)+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(h_{k-\lfloor\frac{k-1}{r}\rfloor
r,i}^{(r-1)}\right)$ (92)
$\displaystyle=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,i+s}^{(r-1)}\right)\,,$

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

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

$\displaystyle\begin{aligned}
h_{k,0}^{(r)}=&\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-2}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,0}^{(r)}\right)+\frac{1}{\left(\lfloor\frac{k-1}{r}\rfloor-1\right)!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor-1}\left(\widehat{h}_{k-\left(\lfloor\frac{k-1}{r}\rfloor-1\right)r,0}^{(r)}\right)\\\
&+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(\widehat{h}_{k-\lfloor\frac{k-1}{r}\rfloor
r,0}^{(r)}\right)\end{aligned}$ (93)

$\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
(88), we have that

$k-jr\geq k-\lfloor\frac{k-1}{r}\rfloor r+2r\geq 2r+1\,,$

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

$k=mr+f\;\;\mathrm{with}\;\;1\leq f\leq
r-1\;\;\mathrm{and}\;\;m\in\mathbb{N}\,.$

From (91) 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 (62) and (63) hold, respectively, for the third and second
term of (93). Thus, we can write equation (93) as

$\displaystyle h_{k,0}^{(r)}=$
$\displaystyle\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-2}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,s}^{(r-1)}\right)+\frac{1}{\left(\lfloor\frac{k-1}{r}\rfloor-1\right)!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor-1}\left(h_{k-\left(\lfloor\frac{k-1}{r}\rfloor-1\right)r,0}^{(r-1)}\right)\,,$

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

  * (ii)

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

$k=mr\;\;\mathrm{with}\;\;m\in\mathbb{N}\,.$

From (90) 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
(62) and (64) hold, respectively, for the third and the second term of (93).
Thus, we can write Eq (93) as

$\displaystyle h_{k,0}^{(r)}=$
$\displaystyle\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,s}^{(r-1)}\right)\,,$ (94)

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

$\displaystyle
h_{k,1}^{(r)}=\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,1}^{(r)}\right)+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(\widehat{h}_{k-\lfloor\frac{k-1}{r}\rfloor
r,1}^{(r)}\right)$ (95)

$\forall\>k\geq r+1\,,k\neq m\,r\,(m\in\mathbb{N})\,$. In view of the
inequalities (88) and (91) we then readily find that the definitions (67) and
(65) hold, respectively, for the first and second part of the previous
equation, i.e., (95) leads to

$\displaystyle h_{k,1}^{(r)}=$
$\displaystyle\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,1+s}^{(r-1)}\right)\,,$ (96)

$\forall\>k\geq r+1\,,k\neq m\,r\,(m\in\mathbb{N})\,$. Finally, recalling
again (88) and (90), we can write Eq (79) as555In the last equality, due to
(90), we have that if $j=\lfloor(k-1)/r\rfloor\,$, then
$\lfloor(k-jr-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.

$\displaystyle\begin{aligned}
h_{k,1}^{(r)}=&\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\frac{1}{j!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}\left(\widehat{h}_{k-jr,1}^{(r)}\right)+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(\widehat{h}_{k-\lfloor\frac{k-1}{r}\rfloor
r,1}^{(r)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)\\\
\phantom{h_{k,1}^{(r)}}=&\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor-1}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,1+s}^{(r-1)}\right)+\frac{1}{\lfloor\frac{k-1}{r}\rfloor!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{\lfloor(k-1)/r\rfloor}\left(h_{k-\lfloor\frac{k-1}{r}\rfloor
r,1}^{(r-1)}\right)\\\
&+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)=\\\
=&\sum_{j=0}^{\lfloor\frac{k-1}{r}\rfloor}\sum_{s=0}^{\lfloor\frac{k-jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{k-jr-
sr,1+s}^{(r-1)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)\,,\end{aligned}$
(97)

$\forall\>k\geq r+1\,,k=m\,r\,\,(m\in\mathbb{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

$\begin{split}&\sum_{j=0}^{m-1}\sum_{s=0}^{\lfloor\frac{mr-
jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{mr-
jr-
sr,1+s}^{(r-1)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)=\\\
&\sum_{j=0}^{m-2}\sum_{s=0}^{\lfloor\frac{mr-
jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{mr-
jr-
sr,1+s}^{(r-1)}\right)+\frac{1}{(m-1)!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(h_{r,1}^{(r-1)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m}\left({\bm{\omega}}\cdot{\bm{p}}\right)=\\\
&\sum_{j=0}^{m-2}\sum_{s=0}^{\lfloor\frac{mr-
jr-1}{r}\rfloor}\frac{1}{j!\,s!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{j}L_{\chi_{1}^{(r)}({\bm{q}})}^{s}\left(h_{mr-
jr-
sr,1+s}^{(r-1)}\right)+\frac{1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(h_{r,1}^{(r-1)}+L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}\left({\bm{\omega}}\cdot{\bm{p}}\right)\right)\\\
&+\frac{m-1}{m!}L_{\chi_{2}^{(r)}({\bm{q}},{\bm{p}})}^{m-1}\left(h_{r,1}^{(r-1)}\right)\,.\end{split}$
(98)

Use of Eq (66) and the homological equation (70) then concludes the proof.∎

##  4 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]).

## References

  * [1] M. Bartuccelli, G. Gentile, Lindstedt series for perturbations of isochronous systems: a review of the general theory, Rev. Math. Phys., 14 (2002), 121–171. http://doi.org/10.1142/S0129055X02001120
  * [2] H. Christodoulidi, C. Efthymiopoulos, T. Bountis, Energy localization on q-tori, long-term stability, and the interpretation of fermi-pasta-ulam recurrences, Phys. Rev. E, 81 (2010), 016210. http://doi.org/10.1103/PhysRevE.81.016210
  * [3] H. Christodoulidi, C. Efthymiopoulos, Low-dimensional $q$-tori in FPU lattices: dynamics and localization properties, Physica D, 261 (2013), 92–113. http://doi.org/10.1016/j.physd.2013.07.007
  * [4] L. Corsi, G. Gentile, M. Procesi, KAM theory in configuration space and cancellations in the Lindstedt series, Commun. Math. Phys., 302 (2011), 359–402. http://doi.org/10.1007/s00220-010-1131-7
  * [5] C. Efthymiopoulos, A. Giorgilli, G. Contopoulos, Nonconvergence of formal integrals II. Improved estimates for the optimal order of truncation, J. Phys. A: Math. Gen., 37 (2004), 10831. http://doi.org/10.1088/0305-4470/37/45/008
  * [6] L. H. Eliasson, Absolutely convergent series expansions for quasi periodic motions, Mathematical Physics Electronic Journal, 2 (1996), 1–33.
  * [7] S. Flach, M. V. Ivanchenko, O. I. Kanakov, q-Breathers and the fermi-pasta-ulam problem, Phys. Rev. Lett., 95 (2005), 064102. http://doi.org/10.1103/PhysRevLett.95.064102
  * [8] S. Flach, M. V. Ivanchenko, O. I. Kanakov, $q$-breathers in Fermi-Pasta-Ulam chains: existence, localization, and stability, Phys. Rev. E, 73 (2006), 036618. http://doi.org/10.1103/PhysRevE.73.036618
  * [9] M. Hénon, C. Heiles, The applicability of the third integral of motion: some numerical experiments, The astronomical journal, 69, 73, 1964.
  * [10] L. Galgani, Foundations of physics in Milan, Padua and Paris. Newtonian trajectories from celestial mechanics to atomic physics, Mathematics in Engineering, 3 (2021), 1–24. http://doi.org/10.3934/mine.2021045
  * [11] G. Gallavotti, Twistless KAM tori, Commun. Math. Phys., 164 (1994), 145–156. http://doi.org/10.1007/BF02108809
  * [12] G. Gallavotti, Twistless KAM tori, quasi flat homoclinic intersections, and other cancellations in the perturbation series of certain completely integrable Hamiltonian systems. A review, Rev. Math. Phys., 6 (1994), 343–411. https://doi.org/10.1142/S0129055X9400016X
  * [13] A. Giorgilli, Small denominators and exponential stability: from Poincaré to the present time, Seminario Mat. e. Fis. di Milano, 68 (1998), 19. http://doi.org/10.1007/BF02925829
  * [14] A. Giorgilli, U. Locatelli, On classical series expansions for quasi-periodic motions, Mathematical Physics Electronic Journal, 3 (1997), 1–25.
  * [15] A. Jorba, J. Masdemont, Dynamics in the center manifold of the collinear points of the restricted three body problem, Physica D, 132 (1999), 189–213. http://doi.org/10.1016/S0167-2789(99)00042-1
  * [16] A. Lindstedt, Ueber die allgemeine form der integrale des dreikörperproblems, Astron. Nachr., 105 (1883), 97–112. http://doi.org/10.1002/ASNA.18831050702
  * [17] A. Lindstedt, Ueber die integration einer gewissen differentialgleichung, Astron. Nachr., 104 (1883), 145–150. http://doi.org/10.1002/asna.18831041002
  * [18] A. Lindstedt, Sur les séries trigonométriques dans le problème des trois corps, Bulletin Astronomique, Serie I, 3 (1886), 217–221.
  * [19] H. Poincaré, Les méthodes nouvelles de la mécanique céleste. Tome II, Gauthier-Villars et fils (Paris) 1893, Méthodes de MM. Newcomb, Gyldén, Lindstedt et Bohlin.
  * [20] M. Sansottera, V. Danesi, Kolmogorov variation: Kam with knobs (à la Kolmogorov), 2021, arXiv:2109.06345.

##  Appendix 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$

$\mathcal{H}(x,p)=\frac{\omega_{0}}{2}\left(p^{2}+x^{2}\right)-\varepsilon\frac{x^{3}}{3}\,.$
(99)

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

$\mathcal{H}(q,J)=\omega_{0}\,J-\frac{\varepsilon}{\sqrt{2}}J^{3/2}\sin(q)+\frac{\varepsilon}{3\,\sqrt{2}}J^{3/2}\sin(3q)\,.$
(100)

Finally, we define the translation $J=J_{0}+p\,$ leading to (apart from
constants)

$\begin{split}\mathcal{H}(q,p)=\omega_{0}\,p-\frac{\varepsilon}{\sqrt{2}}(J_{0}+p)^{3/2}\sin(q)+\frac{\varepsilon}{3\,\sqrt{2}}(J_{0}+p)^{3/2}\sin(3q)\,.\end{split}$
(101)

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 (101) we perform an expansion in $p$
leading to (apart from constants)

$\begin{split}\mathcal{H}(q,p)&=\omega_{0}\,p-\frac{\varepsilon\,J_{0}^{3/2}}{\sqrt{2}}\,\sin(q)-\frac{3\,\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(q)-\frac{3\,\varepsilon}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,p^{2}\,\sin(q)+\frac{\varepsilon}{16\,\sqrt{2}\,J_{0}^{3/2}}\,p^{3}\,\sin(q)\\\
&\phantom{=}+\frac{\varepsilon\,J_{0}^{3/2}}{3\,\sqrt{2}}\,\sin(3q)+\frac{\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(3q)+\frac{\varepsilon}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,p^{2}\,\sin(3q)-\frac{\varepsilon}{48\,\sqrt{2}\,J_{0}^{3/2}}\,p^{3}\,\sin(3q)\,.\end{split}$
(102)

###  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 (100)

$\displaystyle{\dot{q}(t)=\frac{\partial\mathcal{H}}{\partial
J}=\omega_{0}-\frac{3\,\varepsilon\,\sqrt{J}}{2\,\sqrt{2}}\,\sin(q)+\frac{\varepsilon\,\sqrt{J}}{2\,\sqrt{2}}\,\sin(3q)}$
(103a) $\displaystyle{\dot{J}(t)=-\frac{\partial\mathcal{H}}{\partial
q}=\frac{\varepsilon\,J^{3/2}}{\sqrt{2}}\,\cos(q)-\frac{\varepsilon\,J^{3/2}}{\sqrt{2}}\,\cos(3q)}\,.$
(103b)

Replacing, as before, the expressions

$\begin{split}&\omega_{0}=\omega+\varepsilon\,a_{1}+\varepsilon^{2}\,a_{2}+\varepsilon^{3}\,a_{3}\,,\\\
&q(t)=q_{0}(t)+\varepsilon\,q_{1}(t)+\varepsilon^{2}\,q_{2}(t)+\varepsilon^{3}\,q_{3}(t)\,,\\\
&J(t)=J_{0}(t)+\varepsilon
J_{1}(t)+\varepsilon^{2}J_{2}(t)+\varepsilon^{3}J_{3}(t)\,,\end{split}$ (104)

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 (103a) and (103b). At order zero we
have

$\dot{q_{0}}(t)=\omega\,,\qquad\dot{J_{0}}(t)=0\quad\Longrightarrow\quad
q_{0}(t)=\omega\,t\,,\qquad J_{0}(t)=J_{0}\,,$

where we fix the initial phase $q_{0}(0)=0$ and $J_{0}(0)=J_{0}\,$. At order
one, we find

$\begin{split}&\displaystyle{\dot{q_{1}}(t)=a_{1}-\frac{3\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,\sin(\omega\,t)+\frac{\sqrt{J_{0}}}{2\,\sqrt{2}}\,\sin(3\,\omega\,t)},\\\
&\displaystyle{\dot{J_{1}}(t)=\frac{J_{0}^{3/2}}{\sqrt{2}}\,\cos(\omega\,t)-\frac{J_{0}^{3/2}}{\sqrt{2}}\,\cos(3\,\omega\,t)}\,.\end{split}$
(105)

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

$q_{1}(t)=-\frac{2\,\sqrt{2}\,\sqrt{J_{0}}}{3\,\omega}+\frac{3\,\sqrt{J_{0}}}{2\,\sqrt{2}\,\omega}\,\cos(\omega\,t)-\frac{\sqrt{J_{0}}}{6\,\sqrt{2}\,\omega}\,\cos(3\,\omega\,t)\,,\quad
J_{1}(t)=\frac{J_{0}^{3/2}}{\sqrt{2}\,\omega}\,\sin(\omega\,t)-\frac{J_{0}^{3/2}}{3\,\sqrt{2}\,\omega}\,\sin(3\,\omega\,t)$

yielding the constants $q_{1}(0)=J_{1}(0)=0\,$. At order two we now get

$\begin{split}&\displaystyle{\dot{q_{2}}(t)=a_{2}-\frac{5\,J_{0}}{6\,\omega}+\frac{J_{0}}{\omega}\,\cos(\omega\,t)+\frac{3\,J_{0}}{8\,\omega}\,\cos(2\,\omega\,t)-\frac{J_{0}}{\omega}\,\cos(3\,\omega\,t)+\frac{J_{0}}{2\,\omega}\,\cos(4\,\omega\,t)-\frac{J_{0}}{24\,\omega}\,\cos(6\,\omega\,t)},\\\
&\displaystyle{\dot{J_{2}}(t)=\frac{2\,J_{0}^{2}}{3\,\omega}\,\sin(\omega\,t)+\frac{4\,J_{0}^{2}}{3\,\omega}\,\sin(2\,\omega\,t)-\frac{2\,J_{0}^{2}}{\omega}\,\sin(3\,\omega\,t)+\frac{2\,J_{0}^{2}}{3\,\omega}\,\sin(4\,\omega\,t)}\,;\end{split}$
(106)

To compensate for the constant term in $\dot{q_{2}}(t)\,$ we now set
$a_{2}=5\,J_{0}/(6\,\omega)$. Then,

$\begin{split}&\displaystyle{q_{2}(t)=\frac{J_{0}\,}{\omega^{2}}\,\sin(\omega\,t)+\frac{3\,J_{0}}{16\,\omega^{2}}\,\sin(2\,\omega\,t)-\frac{J_{0}}{3\,\omega^{2}}\,\sin(3\,\omega\,t)+\frac{J_{0}}{8\,\omega^{2}}\,\sin(4\,\omega\,t)-\frac{J_{0}}{144\,\omega^{2}}\,\sin(6\,\omega\,t)},\\\
&\displaystyle{J_{2}(t)=\frac{5\,J_{0}^{2}}{6\,\omega^{2}}-\frac{2\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(\omega\,t)-\frac{2\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(2\,\omega\,t)+\frac{2\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(3\,\omega\,t)-\frac{J_{0}^{2}}{6\,\omega^{2}}\,\cos(4\,\omega\,t)}\,.\end{split}$
(107)

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

$\begin{split}&q(t)=\omega\,t-\frac{2\,\sqrt{2}\,\varepsilon\,\sqrt{J_{0}}}{3\,\omega}+\frac{3\,\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}\,\omega}\,\cos(\omega\,t)-\frac{\varepsilon\,\sqrt{J_{0}}}{6\,\sqrt{2}\,\omega}\,\cos(3\,\omega\,t)+\frac{\,\varepsilon^{2}\,J_{0}}{\omega^{2}}\,\sin(\omega\,t)+\frac{3\,\varepsilon^{2}\,J_{0}}{16\,\omega^{2}}\,\sin(2\,\omega\,t)\\\
&\phantom{q(t)=}-\frac{\varepsilon^{2}\,J_{0}}{3\,\omega^{2}}\,\sin(3\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}}{8\,\omega^{2}}\,\sin(4\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}}{144\,\omega^{2}}\,\sin(6\,\omega\,t)-\frac{38\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{81\omega^{3}}+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{\sqrt{2}\,\omega^{3}}\,\cos(\omega\,t)\\\
&\phantom{q(t)=}-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{2\,\sqrt{2}\,\omega^{3}}\,\cos(2\,\omega\,t)+\frac{145\,\varepsilon^{3}\,J_{0}^{3/2}}{144\,\sqrt{2}\,\omega^{3}}\,\cos(3\,\omega\,t)-\frac{\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{3\,\omega^{3}}\,\cos(4\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{16\,\sqrt{2}\,\omega^{3}}\,\cos(5\,\omega\,t)\\\
&\phantom{q(t)=}+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{18\,\sqrt{2}\,\omega^{3}}\,\cos(6\,\omega\,t)-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{48\,\sqrt{2}\,\omega^{3}}\,\cos(7\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{1296\,\sqrt{2}\,\omega^{3}}\,\cos(9\,\omega\,t)\,,\\\\[5.59721pt]
&J(t)=J_{0}+\frac{\varepsilon\,J_{0}^{3/2}}{\sqrt{2}\,\omega}\,\sin(\omega\,t)-\frac{\varepsilon\,J_{0}^{3/2}}{3\,\sqrt{2}\,\omega}\,\sin(3\,\omega\,t)+\frac{5\,\varepsilon^{2}\,J_{0}^{2}}{6\,\omega^{2}}-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(\omega\,t)-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(2\,\omega\,t)\\\
&\phantom{J(t)=}+\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(3\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}^{2}}{6\,\omega^{2}}\,\cos(4\,\omega\,t)+\frac{7\,\varepsilon^{3}\,J_{0}^{5/2}}{4\,\sqrt{2}\,\omega^{3}}\,\sin(\omega\,t)-\frac{8\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega^{3}}\,\sin(2\,\omega\,t)\\\
&\phantom{J(t)=}+\frac{13\,\varepsilon^{3}\,J_{0}^{5/2}}{8\,\sqrt{2}\,\omega^{3}}\,\sin(3\,\omega\,t)-\frac{4\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega^{3}}\,\sin(4\,\omega\,t)+\frac{7\,\varepsilon^{3}\,J_{0}^{5/2}}{72\,\sqrt{2}\,\omega^{3}}\,\sin(5\,\omega\,t)\,,\\\\[5.59721pt]
&\omega=\omega_{0}-\frac{5\,\varepsilon^{2}J_{0}}{6\,\omega}\,.\end{split}$
(108)

###  A.2 Kolmogorov normal form

Starting from the Hamiltonian (102), 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 (102), we have

$\begin{split}\mathcal{H}(q,p)&=\omega\,p+\varepsilon\,a_{1}\,p-\frac{\varepsilon\,J_{0}^{3/2}}{\sqrt{2}}\,\sin(q)-\frac{3\,\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(q)-\frac{3\,\varepsilon}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,p^{2}\,\sin(q)\\\
&\phantom{=}+\frac{\varepsilon}{16\,\sqrt{2}\,J_{0}^{3/2}}\,p^{3}\,\sin(q)+\frac{\varepsilon\,J_{0}^{3/2}}{3\,\sqrt{2}}\,\sin(3q)+\frac{\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(3q)+\frac{\varepsilon}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,p^{2}\,\sin(3q)\\\
&\phantom{=}-\frac{\varepsilon}{48\,\sqrt{2}\,J_{0}^{3/2}}\,p^{3}\,\sin(3q)+\varepsilon^{2}\,a_{2}\,p+\varepsilon^{3}\,a_{3}\,p\,.\end{split}$
(109)

At first order, we have

$\left\langle h_{1,1}^{(0)}\right\rangle_{{q}}=\left\langle
a_{1}\,p-\frac{3\,\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(q)+\frac{\sqrt{J_{0}}}{2\,\sqrt{2}}\,p\,\sin(3\,q)\right\rangle_{{q}}=0\qquad\Longrightarrow\qquad
a_{1}=0\,;$ (110)

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

$h_{1,0}^{(0)}(q)=-\frac{J_{0}^{3/2}}{\sqrt{2}}\sin(q)+\frac{J_{0}^{3/2}}{3\,\sqrt{2}}\,\sin(3\,q)\,,$

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

$X^{(1)}=\frac{J_{0}^{3/2}}{\sqrt{2}\,\omega}\,\cos(q)-\frac{J_{0}^{3/2}}{9\,\sqrt{2}\,\omega}\,\cos(3\,q)\,.$

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

$\widehat{\mathcal{H}}^{(1)}=\exp\left(L_{\varepsilon\,\chi_{1}^{(1)}(q)}\right)\mathcal{H}\,.$

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

$\widehat{h}_{1,1}^{(1)}=-\frac{3\,\sqrt{J_{0}}\,p}{2\,\sqrt{2}}\,\sin(q)+\frac{\sqrt{J_{0}}\,p}{2\,\sqrt{2}}\,\sin(3\,q)\,.$

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

$S^{(1)}=-\frac{2\,\sqrt{2}\,\sqrt{J_{0}}}{3\,\omega}$

and

$\chi_{2}^{(1)}(q,\,p)=-\frac{2\,\sqrt{2}\,\sqrt{J_{0}}\,p}{3\,\omega}+\frac{3\,\sqrt{J_{0}}\,p}{2\,\sqrt{2}\,\omega}\,\cos(q)-\frac{\sqrt{J_{0}}\,p}{6\,\sqrt{2}\,\omega}\,\cos(3\,q)\,.$

Thus, the Hamiltonian

$\mathcal{H}^{(1)}=\exp\left(L_{\varepsilon\,\chi_{2}^{(1)}(q,\,p)}\right)\widehat{\mathcal{H}}^{(1)}\,$

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
$-\displaystyle{\frac{3\,p^{2}}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,\sin(q)}+\displaystyle{\frac{p^{2}}{8\,\sqrt{2}\,\sqrt{J_{0}}}\,\sin(3\,q)}+\displaystyle{\frac{p^{3}}{16\,\sqrt{2}\,J_{0}^{3/2}}\,\sin(q)}-\displaystyle{\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:

$\begin{split}&a_{2}=\frac{5\,J_{0}}{6\,\omega}\,,\\\
&\chi_{1}^{(2)}=-\frac{5\,J_{0}^{2}\,q}{12\,\omega^{2}}+\frac{5\,J_{0}^{2}}{16\,\omega^{2}}\,\sin(2\,q)-\frac{J_{0}^{2}}{16\,\omega^{2}}\,\sin(4\,q)+\frac{J_{0}^{2}}{144\,\omega^{2}}\,\sin(6\,q)\,,\\\
&\chi_{2}^{(2)}=\frac{J_{0}\,p}{2\,\omega^{2}}\,\sin(q)+\frac{13\,J_{0}\,p}{32\,\omega^{2}}\,\sin(2\,q)-\frac{J_{0}\,p}{6\,\omega^{2}}\,\sin(3\,q)+\frac{J_{0}\,p}{288\,\omega^{2}}\,\sin(6\,q)\end{split}$
(111)

and

$\displaystyle\begin{split}&a_{3}=0\,,\\\
&\chi_{1}^{(3)}=\frac{49\,J_{0}^{5/2}}{96\,\sqrt{2}\,\omega^{3}}\,\cos(q)-\frac{5\,J_{0}^{5/2}}{6\,\sqrt{2}\,\omega^{3}}\,\cos(2\,q)+\frac{43\,J_{0}^{5/2}}{432\,\sqrt{2}\,\omega^{3}}\,\cos(3\,q)+\frac{J_{0}^{5/2}}{3\,\sqrt{2}\,\omega^{3}}\,\cos(4\,q)\\\
&\phantom{\chi_{1}^{(3)}=}-\frac{29\,J_{0}^{5/2}}{240\,\sqrt{2}\,\omega^{3}}\,\cos(5\,q)-\frac{J_{0}^{5/2}}{18\,\sqrt{2}\,\omega^{3}}\,\cos(6\,q)+\frac{7\,J_{0}^{5/2}}{192\,\sqrt{2}\,\omega^{3}}\,\cos(7\,q)-\frac{11\,J_{0}^{5/2}}{5184\,\sqrt{2}\,\omega^{3}}\,\cos(9\,q)\,,\\\
&\chi_{2}^{(3)}=-\frac{107\,J_{0}^{3/2}\,p}{162\,\sqrt{2}\,\omega^{3}}+\frac{295\,J_{0}^{3/2}\,p}{192\,\sqrt{2}\,\omega^{3}}\,\cos(q)-\frac{47\,J_{0}^{3/2}\,p}{36\,\sqrt{2}\,\omega^{3}}\,\cos(2\,q)+\frac{133\,J_{0}^{3/2}\,p}{288\,\sqrt{2}\,\omega^{3}}\,\cos(3\,q)\\\
&\phantom{\chi_{2}^{(3)}=}-\frac{J_{0}^{3/2}\,p}{18\,\sqrt{2}\,\omega^{3}}\,\cos(4\,q)+\frac{13\,J_{0}^{3/2}\,p}{288\,\sqrt{2}\,\omega^{3}}\,\cos(5\,q)-\frac{J_{0}^{3/2}\,p}{36\,\sqrt{2}\,\omega^{3}}\,\cos(6\,q)+\frac{7\,J_{0}^{3/2}\,p}{1152\,\sqrt{2}\,\omega^{3}}\,\cos(7\,q)\\\
&\phantom{\chi_{2}^{(3)}=}-\frac{J_{0}^{3/2}\,p}{10368\,\sqrt{2}\,\omega^{3}}\,\cos(9\,q)\,,\end{split}$
(112)

for the generating functions. The final Hamiltonian is

$\mathcal{H}^{(3)}(\tilde{q},\tilde{p})=\omega\,\tilde{p}+\varepsilon\,R_{1}(q,\,p)+\varepsilon^{2}\,R_{2}(q,\,p)+\varepsilon^{3}\,R_{3}(q,\,p)$

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

$\displaystyle
q(t)=\omega\,t-\frac{2\,\sqrt{2}\,\varepsilon\,\sqrt{J_{0}}}{3\,\omega}+\frac{3\,\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}\,\omega}\,\cos(\omega\,t)-\frac{\varepsilon\,\sqrt{J_{0}}}{6\,\sqrt{2}\,\omega}\,\cos(3\,\omega\,t)+\frac{\,\varepsilon^{2}\,J_{0}}{\omega^{2}}\,\sin(\omega\,t)+\frac{3\,\varepsilon^{2}\,J_{0}}{16\,\omega^{2}}\,\sin(2\,\omega\,t)$
$\displaystyle\phantom{q(t)=}-\frac{\varepsilon^{2}\,J_{0}}{3\,\omega^{2}}\,\sin(3\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}}{8\,\omega^{2}}\,\sin(4\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}}{144\,\omega^{2}}\,\sin(6\,\omega\,t)-\frac{38\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{81\omega^{3}}+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{\sqrt{2}\,\omega^{3}}\,\cos(\omega\,t)$
$\displaystyle\phantom{q(t)=}-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{2\,\sqrt{2}\,\omega^{3}}\,\cos(2\,\omega\,t)+\frac{145\,\varepsilon^{3}\,J_{0}^{3/2}}{144\,\sqrt{2}\,\omega^{3}}\,\cos(3\,\omega\,t)-\frac{\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{3\,\omega^{3}}\,\cos(4\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{16\,\sqrt{2}\,\omega^{3}}\,\cos(5\,\omega\,t)$
$\displaystyle\phantom{q(t)=}+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{18\,\sqrt{2}\,\omega^{3}}\,\cos(6\,\omega\,t)-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{48\,\sqrt{2}\,\omega^{3}}\,\cos(7\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{1296\,\sqrt{2}\,\omega^{3}}\,\cos(9\,\omega\,t)\,,$
$\displaystyle
p(t)=\frac{\varepsilon\,J_{0}^{3/2}}{\sqrt{2}\,\omega}\,\sin(\omega\,t)-\frac{\varepsilon\,J_{0}^{3/2}}{3\,\sqrt{2}\,\omega}\,\sin(3\,\omega\,t)+\frac{5\,\varepsilon^{2}\,J_{0}^{2}}{6\,\omega^{2}}-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(\omega\,t)-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(2\,\omega\,t)$
$\displaystyle\phantom{J(t)=}+\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega^{2}}\,\cos(3\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}^{2}}{6\,\omega^{2}}\,\cos(4\,\omega\,t)+\frac{7\,\varepsilon^{3}\,J_{0}^{5/2}}{4\,\sqrt{2}\,\omega^{3}}\,\sin(\omega\,t)-\frac{8\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega^{3}}\,\sin(2\,\omega\,t)$
$\displaystyle\phantom{J(t)=}+\frac{13\,\varepsilon^{3}\,J_{0}^{5/2}}{8\,\sqrt{2}\,\omega^{3}}\,\sin(3\,\omega\,t)-\frac{4\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega^{3}}\,\sin(4\,\omega\,t)+\frac{7\,\varepsilon^{3}\,J_{0}^{5/2}}{72\,\sqrt{2}\,\omega^{3}}\,\sin(5\,\omega\,t)\,,$
$\displaystyle J(t)=J_{0}+p(t)\,.$ (113)

Recalling also the computed values of $a_{1}\,$, $a_{2}\,$ and $a_{3}$ (Eqs
(110)–(112)) we have also

$\omega=\omega_{0}-\frac{5\,\varepsilon^{2}J_{0}}{6\,\omega}\,.$

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

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:

$\displaystyle
q(t)=\omega\,t-\frac{2\,\sqrt{2}\,\varepsilon\,\sqrt{J_{0}}}{3\omega_{0}}+\frac{3\,\varepsilon\,\sqrt{J_{0}}}{2\,\sqrt{2}\,\omega_{0}}\,\cos(\omega\,t)-\frac{\varepsilon\,\sqrt{J_{0}}}{6\,\sqrt{2}\,\omega_{0}}\,\cos(3\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}}{\omega_{0}^{2}}\,\sin(\omega\,t)+\frac{3\,\varepsilon^{2}\,J_{0}}{16\,\omega_{0}^{2}}\,\sin(2\,\omega\,t)$
$\displaystyle\phantom{q(t)=}-\frac{\varepsilon^{2}\,J_{0}}{3\,\omega_{0}^{2}}\,\sin(3\,\omega\,t)+\frac{\varepsilon^{2}\,J_{0}}{8\,\omega_{0}^{2}}\,\sin(4\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}}{144\,\omega_{0}^{2}}\,\sin(6\,\omega\,t)-\frac{83\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{81\,\omega_{0}^{3}}+\frac{9\,\varepsilon^{3}\,J_{0}^{3/2}}{4\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(\omega\,t)$
$\displaystyle\phantom{q(t)=}-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{2\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(2\,\omega\,t)+\frac{125\,\varepsilon^{3}\,J_{0}^{3/2}}{144\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(3\omega\,t)-\frac{\sqrt{2}\,\varepsilon^{3}\,J_{0}^{3/2}}{3\,\omega_{0}^{3}}\,\cos(4\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{16\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(5\,\omega\,t)$
$\displaystyle\phantom{q(t)=}+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{18\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(6\,\omega\,t)-\frac{\varepsilon^{3}\,J_{0}^{3/2}}{48\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(7\,\omega\,t)+\frac{\varepsilon^{3}\,J_{0}^{3/2}}{1296\,\sqrt{2}\,\omega_{0}^{3}}\,\cos(9\,\omega\,t)\,,$
$\displaystyle
p(t)=\frac{\varepsilon\,J_{0}^{3/2}}{\sqrt{2}\,\omega_{0}}\,\sin(\omega\,t)-\frac{\varepsilon\,J_{0}^{3/2}}{3\,\sqrt{2}\,\omega_{0}}\,\sin(3\,\omega\,t)+\frac{5\,\varepsilon^{2}\,J_{0}^{2}}{6\,\omega_{0}^{2}}-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega_{0}^{2}}\,\cos(\omega\,t)-\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega_{0}^{2}}\,\cos(2\,\omega\,t)$
$\displaystyle\phantom{p(t)=}+\frac{2\,\varepsilon^{2}\,J_{0}^{2}}{3\,\omega_{0}^{2}}\,\cos(3\,\omega\,t)-\frac{\varepsilon^{2}\,J_{0}^{2}}{6\,\omega_{0}^{2}}\,\cos(4\,\omega\,t)+\frac{31\,\varepsilon^{3}\,J_{0}^{5/2}}{12\,\sqrt{2}\,\omega_{0}^{3}}\,\sin(\omega\,t)-\frac{8\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega_{0}^{3}}\,\sin(2\,\omega\,t)$
$\displaystyle\phantom{p(t)=}+\frac{97\,\varepsilon^{3}\,J_{0}^{5/2}}{72\,\sqrt{2}\,\omega_{0}^{3}}\,\sin(3\,\omega\,t)-\frac{4\,\sqrt{2}\,\varepsilon^{3}\,J_{0}^{5/2}}{9\,\omega_{0}^{3}}\,\sin(4\,\omega\,t)+\frac{7\,\varepsilon^{3}\,J_{0}^{5/2}}{72\,\sqrt{2}\,\omega_{0}^{3}}\,\sin(5\,\omega\,t)\,,$
$\displaystyle\omega=\omega_{0}-\frac{5\,\varepsilon^{2}\,J_{0}}{6\,\omega_{0}}\,.$
(114)

##  Appendix B Proof of Lemma 1

Consider the Hamiltonian:

$\mathcal{H}(q,J)=\omega_{0}\,J+\varepsilon\,h(q,J)\,,$

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

$\displaystyle
f(p)=\varepsilon\,(J_{0}+p)^{\frac{k}{2}}=\varepsilon\,J_{0}^{\frac{k}{2}}+\varepsilon\,\frac{k}{2}\,J_{0}^{\frac{k}{2}-1}p+\ldots+\frac{\varepsilon}{n!}\,\frac{k}{2}\left(\frac{k}{2}-1\right)\cdots\left(\frac{k}{2}-n+1\right)J_{0}^{\frac{k}{2}-n}p^{n}\,.$

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

$\displaystyle
L_{\varepsilon\chi_{1}^{(1)}}^{n-1}\left(\varepsilon\,J_{0}^{\frac{k}{2}-n}p^{n}\right)$
$\displaystyle=L_{\varepsilon\chi_{1}^{(1)}}^{n-2}\underbrace{\\{\varepsilon\,J_{0}^{\frac{k}{2}-n}p^{n},\varepsilon\chi_{1}^{(1)}(q)\\}}_{\sim\,\varepsilon^{2}J_{0}^{k-n}\,p^{n-1}}=L_{\varepsilon\chi_{1}^{(1)}}^{n-3}\underbrace{\\{\\{\varepsilon\,J_{0}^{\frac{k}{2}-n}p^{n},\varepsilon\chi_{1}^{(1)}(q)\\},\varepsilon\chi_{1}^{(1)}(q)\\}}_{\sim\,\varepsilon^{3}J_{0}^{\frac{3}{2}\,k-n}\,p^{n-2}}\simeq\underbrace{\ldots\ldots}_{(n-3)\mathrm{times}}$
(115)
$\displaystyle\simeq\,\varepsilon^{n}J_{0}^{\frac{k}{2}-n}\left(J_{0}^{\frac{k}{2}}\right)^{n-1}p=\varepsilon^{n}J_{0}^{\frac{(k-2)n}{2}}\,p\,,$

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
(115) 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

$\displaystyle
L_{\varepsilon\chi_{2}^{(1)}}^{n-1}\left(\varepsilon\,J_{0}^{\frac{k}{2}-1}p\right)$
$\displaystyle=L_{\varepsilon\chi_{2}^{(1)}}^{n-2}\underbrace{\\{\varepsilon\,J_{0}^{\frac{k}{2}-1}p,\varepsilon\chi_{2}^{(1)}(q,p)\\}}_{\sim\,\varepsilon^{2}J_{0}^{k-2}\,p}=\underbrace{\ldots\ldots}_{(n-2)\mathrm{times}}\simeq\,\varepsilon^{n}J_{0}^{\frac{k}{2}-1}\left(J_{0}^{\frac{k}{2}-1}\right)^{n-1}p=\varepsilon^{n}J_{0}^{\frac{(k-2)n}{2}}\,p\,,$

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

$\displaystyle
L_{\varepsilon\chi_{1}^{(1)}}^{n-1}\left(\varepsilon\,J_{0}^{\frac{k}{2}-n+1}p^{n-1}\right)$
$\displaystyle=L_{\varepsilon\chi_{1}^{(1)}}^{n-2}\underbrace{\\{\varepsilon\,J_{0}^{\frac{k}{2}-n+1}p^{n-1},\varepsilon\chi_{1}^{(1)}(q)\\}}_{\sim\,\varepsilon^{2}J_{0}^{k-n+1}\,p^{n-2}}=\underbrace{\ldots\ldots}_{(n-2)\mathrm{times}}\simeq\,\varepsilon^{n}J_{0}^{\frac{k}{2}-n+1}\left(J_{0}^{\frac{k}{2}}\right)^{n-1}p$
$\displaystyle=\varepsilon^{n}J_{0}^{\frac{(k-2)n+2}{2}}\,p\,,$

implying $\chi_{1}^{(n)}\sim\,J_{0}^{\frac{(k-2)n+2}{2}}\,$. This concludes
the proof of the lemma. ∎