life of the ISS is connected to the speed of this diffusion, which eventually
drives the system to the instability. It should be stressed that, besides the
speed of diffusion, the lifetime of the inner orbits also depends on the
initial distance of the system from the instability boundary defined by the
resonance $g_{1}-g_{5}$. This geometric aspect includes the stabilizing role
of general relativity [5, 6], which moves the system away from the instability
boundary by 0.43${}^{\prime\prime}\,{\textrm{yr}}^{-1}$, and the destabilizing
effect of terms of degree 6 in eccentricities and inclinations of the planets
[7].

##  VII Discussion

This work introduces a framework that naturally justifies the statistical
stability shown by the ISS over a timescale comparable to its age. Considering
a forced secular model of the inner planet orbits, the computation of the
Lyapunov spectrum indicates the existence of very different dynamical
timescales. Using the computer algebra system TRIP, we systematically analyze
the Fourier harmonics of the Hamiltonian that become resonant along a
numerically integrated orbital solution spanning 5 Gyr. We uncover three
symmetries that characterize the strongest resonances and that are broken by
weak resonant interactions. These quasi-symmetries generate three QIs of
motion that represent slow variables of the secular dynamics. The size of the
leading symmetry-breaking resonances suggests that the QIs are related to the
smallest Lyapunov exponents. The claim that the QIs are among the slowest
d.o.f. of the dynamics constitutes the central point of this work. On the one
hand, it is supported by the analysis of the underlying Hamiltonian
$\mathbb{H}_{2n}^{\bullet}$, in which one neglects the forcing mainly exerted
by Saturn, Uranus, and Neptune and, as a consequence, the diffusion of the QIs
is greatly reduced. On the other hand, the geometric framework established by
the PCA of the orbital solutions independently confirms that the QIs are
statistically the slowest linear variables of the dynamics. We give strong
evidence that the QIs of motion play a critical role in the statistical
stability of the ISS over the Solar System lifetime, by adiabatically
constraining the long-term chaotic diffusion of the orbits.

###  VII.1 Inner Solar System among classical quasi-integrable systems

It is valuable to contextualize the dynamics of the ISS in the class of
classical quasi-integrable systems. A comparison with the Fermi-Pasta-Ulam-
Tsingou problem, in particular, deserves to be made. This concerns the
dynamics of a one-dimensional chain of identical masses coupled by nonlinear
springs. For weak nonlinearity, the normal modes of oscillation remain far
from the energy equipartition expected from statistical mechanics for a very
long time [13]. One way to explain the lack of energy equipartition reported
by Fermi and collaborators is through the closeness of the FPUT problem to the
integrable Toda dynamics [63, 64, 65]. This translates in a very slow
thermalization of the action variables of the Toda problem and of the
corresponding integrals of motion along the FPUT flow [65, 66, 67, 15, 68,
69]. In the framework of the present study, the very long dynamical half-life
of the ISS is also likely to be the result of the slow diffusion of some
dynamical quantities, the QIs of motion. We find, in particular, an underlying
Hamiltonian $\mathbb{H}_{2n}^{\bullet}$ for which this diffusion is greatly
reduced, as a consequence of neglecting the forcing mainly exerted by Saturn,
Uranus, and Neptune. This results in a dynamics that can be considered as
stable in an astronomical sense. We stress that, differently from the FPUT
problem, $\mathbb{H}_{2n}^{\bullet}$ is not integrable as the Toda
Hamiltonian. It is indeed chaotic and shares with the original Hamiltonian
$\mathbb{H}_{2n}$ the leading Lyapunov exponents. The QIs that we find in this
work are only a small number of functions of the action-angle variables of the
integrable LL dynamics, and are related to the smallest Lyapunov exponents of
the dynamics. Our study suggests that in the FPUT problem the very slow
thermalization occurring beyond the Lyapunov time might be understood in terms
of combinations of the Toda integrals of motion diffusing over very different
timescales.

The long-term diffusion in chaotic quasi-integrable systems should be
generally characterized by a broad range of timescales that results from the
progressive, hierarchical breaking of the symmetries of the underlying
integrable problem by resonant interactions [70, 71, 72]. A hierarchy of
Lyapunov exponents spanning several orders of magnitude, in particular, should
be common among this class of systems [e.g., 73].

###  VII.2 Methods

The long-term dynamics of the ISS is described by a moderate but not small
number of d.o.f., which places it far from the typical application fields of
celestial mechanics and statistical physics. The first discipline often
studies dynamical models with very few degrees of freedom, while the second
one deals with the limit of a very large number of bodies. Chaos also requires
a statistical description of the inner planet orbits. But the lack of a
statistical equilibrium, resulting from a slow but ceaseless diffusion of the
system, places the ISS outside the standard framework of ergodic theory. The
kind of approach we develop in this work is heavily based on computer algebra,
in terms of systematic series expansion of the Hamiltonian, manipulation of
the truncated equations of motion, extraction of given Fourier harmonics,
retrieval of polynomial roots, etc. [4, 19]. This allows us to introduce QIs
of motion in a 16-dimensional dynamics by analyzing how action-space
symmetries are progressively broken by resonant interactions. Our effective
method based on the time statistics of resonances arising along a single, very
long numerical integration is alternative to formal approaches that define QIs
via series expansions [e.g. 74, 75]. The practical usefulness of these formal
expansions for a dynamics that covers an intricate, high-dimensional network
of resonances seems indeed doubtful. Through the retrieval of the half-widths
of the symmetry-breaking resonances, computer algebra also permits us to
extend the correspondence between the Lyapunov spectrum and the spectrum of
resonances well beyond the standard relation linking the Lyapunov time to the
strongest resonances888In this regard, it should be noted that a relation
between QIs and Lyapunov exponents has already been highlighted in simple
systems [83, 43].

In the context of dynamical systems with a number of d.o.f. that is not small,
this work also considers an approach based on PCA. The role of this
statistical technique can be twofold. We use PCA as an independent test to
systematically validate the slowness of the QIs. While being introduced semi-
analytically as dynamical quantities that are not affected by the leading
resonances, they can indeed be related to the last PCs. By extension, the
first PCs should probe the directions of the main resonances. This leads to a
second potential application of the PCA, which should offer a way to retrieve
the principal resonant structure of a dynamical system. In this sense, PCA
represents a tool to systematically probe numerical integrations of a complex
dynamics and distill important hidden insights. We emphasize that PCA is the
most basic linear technique of dimensionality reduction and belongs to the
more general class of the unsupervised learning algorithms. There are more
sophisticated methods of feature extraction that can be more robust [e.g. 77,
78] and can incorporate nonlinearity [79]. These methods are often less
intuitive to understand, less straightforward to apply, and harder to
interpret than PCA. Yet, they might be more effective and worth pursuing for
future works.

With long-term numerical integration and a computer algebra system at one’s
disposal, the entire strategy we develop in this work can in principle be
applied to other planetary systems and quasi-integrable Hamiltonian dynamics
with a moderate number of d.o.f.

###### Acknowledgements.

The authors thank M. Gastineau for his assistance with TRIP. F.M. is supported
by a grant of the French Agence Nationale de la Recherche (AstroMeso
ANR-19-CE31-0002-01). N.H.H. is supported by a Ph.D. scholarship of the
Capital Fund Management (CFM) Foundation for Research. This project has been
supported by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (Advanced Grant AstroGeo-885250).
This work was granted access to the HPC resources of MesoPSL financed by the
Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01)
of the program Investissements d’Avenir supervised by the Agence Nationale
pour la Recherche.

##  Appendix A Lyapunov spectrum

(a)  $\mathcal{H}_{4}$: Single orbital solution (b)  $\mathcal{H}_{4}$:
Relative numerical errors

ciao

Figure 8: (a) Positive FT-LCEs of Hamiltonian $\mathcal{H}_{4}$ and
corresponding characteristic timescales for a single initial condition and an
ensemble of 150 random sets of initial tangent vectors. The bands represent
the [5th, 95th] percentile range of the marginal PDFs. The lines denote the
distribution medians. (b) Medians of the relative numerical errors
$\epsilon_{i}$ on the FT-LCEs $\lambda_{i}$, as defined in Eq. (33), for the
ensemble of 150 orbital solution of Fig. 1a.

##### Convergence.

We perform two tests to address the convergence of our implementation of the
Benettin  _et al._ [43] method. We first compute the FT-LCEs for a single
initial condition of $\mathcal{H}_{4}$ and an ensemble of 150 different random
sets of initial tangent vectors. Figure 7a shows the [5th, 95th] percentile
range of the resulting marginal distributions of the positive FT-LCEs over a
time span of 10 Gyr. The distributions shrink with increasing time, eventually
collapsing on single time-dependent values. In this asymptotic regime, the
Benettin  _et al._ [43] algorithm loses memory of the initial tangent vectors
and purely retrieves the FT-LCEs as defined in Eq. (10). Therefore, Fig. 1a
shows asymptotically the dependence of the FT-LCEs on the initial condition
$\bm{z}_{0}$ and represents their statistical distribution over the phase-
space domain explored by the dynamics in a non-ergodic way. The convergence of
the computation is clearly slower for smaller exponents, but a comparison with
Fig. 1a indicates that, even in the case of $\lambda_{8}$, the numerical
uncertainty on the FT-LCEs of each orbital solution at 10 Gyr is negligible
with respect to the width of their ensemble distributions.

To quantitatively estimate the numerical precision on the computed FT-LCEs, we
exploit the symmetry of the spectrum stated in Eq. (9). For a single orbital
solution, the relative numerical error on each exponent $\lambda_{i}$ can be
estimated as

$\epsilon_{i}=\left|\frac{\Delta\lambda_{i}}{\lambda_{i}}\right|.$ (33)

We plot in Fig. 7b the medians of $\epsilon_{i}$ for the ensemble of 150
orbital solutions of Fig. 1a. The relative errors decrease asymptotically with
time, as expected. Even in the case of the smallest exponent, $\lambda_{8}$,
the median error is less than 10% at 10 Gyr.

##### Hamiltonian $\mathcal{H}_{6}$.

We compute for comparison the FT-LCEs of the forced ISS truncated at degree 6
in eccentricities and inclinations, that is, $\mathcal{H}_{6}$. We consider
150 stable orbital solutions with initial conditions very close to the nominal
values of Gauss’s dynamics and random sets of initial tangent vectors, as we
do for the truncation at degree 4. Figure 9 shows the [5th, 95th] percentile
range of the marginal PDF of each FT-LCE estimated from the ensemble of
solutions. Apart from being somewhat larger, the asymptotic distributions of
the exponents are very similar to those of $\mathcal{H}_{4}$ shown in Fig. 1a.

##  Appendix B Vectors $\bgamma_{1}$, $\bgamma_{2}$, $\bgamma_{3}$

We report here the explicit expressions of the vectors
$(\bm{\gamma}_{i})_{i=1}^{3}$. We first give the components of the vector
$\bm{\omega}_{\scaleto{\textsc{LL}}{4pt}}$ of the fundamental precession
frequencies of the inner orbits in the forced Laplace-Lagrange dynamics
(including the leading correction of general relativity) [4]:

$\displaystyle\bm{\omega}_{\scaleto{\textsc{LL}}{4pt}}$
$\displaystyle=(\bm{g}_{\scaleto{\textsc{LL}}{4pt}},\bm{s}_{\scaleto{\textsc{LL}}{4pt}})\approx$
(34) $\displaystyle(5.87,7.46,17.4,18.1,-5.21,-6.59,-18.8,-17.7),$

in units of ${\textrm{arcsec yr}}^{-1}$ (see [80, 81, 82] for comparison with
the frequencies of the Laplace-Lagrange dynamics of the entire Solar System).
One then has

$\displaystyle\bm{\gamma}_{1}$
$\displaystyle=(\bm{0}_{4},\bm{1}_{4})=(0,0,0,0,1,1,1,1),$ (35)
$\displaystyle\bm{\gamma}_{2}$ $\displaystyle=(0,0,-1,-1,1,1,2,2),$
$\displaystyle\bm{\gamma}_{3}$
$\displaystyle=-\bm{\omega}_{\scaleto{\textsc{LL}}{4pt}}+g_{5}\bm{1}_{8}\approx$
$\displaystyle(-1.61,-3.20,-13.2,-13.9,9.47,10.8,23.0,22.0),$

with the components of $\bm{\gamma}_{3}$ in units of ${\textrm{arcsec
yr}}^{-1}$. We recall that $g_{5}\approx
4.257$${}^{\prime\prime}\,{\textrm{yr}}^{-1}$ is a constant in the forced
model of the ISS. The corresponding unit vectors
$(\widehat{\bm{\gamma}}_{i})_{i=1}^{3}$ are given by

$\displaystyle\widehat{\bm{\gamma}}_{1}$ $\displaystyle=(0,0,0,0,1,1,1,1)/2,$
(36) $\displaystyle\widehat{\bm{\gamma}}_{2}$
$\displaystyle=(0,0,-1,-1,1,1,2,2)/2\sqrt{3},$
$\displaystyle\widehat{\bm{\gamma}}_{3}$
$\displaystyle\approx(-0.04,-0.08,-0.33,-0.35,0.24,0.27,0.58,0.55).$

Since $1/2\sqrt{3}\approx 0.289$, the components of
$\widehat{\bm{\gamma}}_{3}$ are only a few percent away from those of
$\widehat{\bm{\gamma}}_{2}$. Therefore, along stable orbital solutions with
typical bounded variations of the Mercury-dominated action variable
$\mathfrak{X}_{1}$, the two quantities $\mathcal{C}_{2}$ and
$\mathcal{E}_{2n}$ exhibit very similar time evolutions. This is not the case
anymore when Mercury orbit reaches high eccentricities.

Figure 9: Positive FT-LCEs $\lambda_{i}$ of Hamiltonian $\mathcal{H}_{6}$ and
corresponding characteristic timescales $\lambda_{i}^{-1}$. The bands
represent the [5th, 95th] percentile range of the marginal PDFs estimated from
an ensemble of 150 stable orbital solutions with very close initial
conditions. The lines denote the distribution medians.

##  Appendix C Ensemble distributions of the quasi-integrals over time

To retrieve the long-term statistical behavior of the QIs, we consider the
ensembles of 1080 numerical integrations of the dynamical models
$\mathcal{H}_{4}$ and $\mathcal{H}_{6}$, with very close initial conditions
and spanning 100 Gyr in the future, that have been presented in Ref. [7]. We
also consider the similar ensembles of solutions for the simplified
Hamiltonians $\mathbb{H}_{4}^{\bullet}$ and $\mathbb{H}_{6}^{\bullet}$ that we
introduce in Sect. IV.4. We report in Fig. 10 the time evolution of the
ensemble PDFs of the low-pass filtered dimensionless QIs and dimensionless
actions $\mathfrak{X}_{1},\Psi_{3}$ for the different models (the cutoff
frequency of the time filter is set to 1 Myr-1, as in Sec. IV.2). More
precisely, to highlight the growth of the statistical dispersion, we consider
at each time the PDF of the signed deviation from the ensemble mean, so that
all the plotted distributions have a null mean. At each time, the PDF
estimation takes into account only the stable orbital solutions, that is,
those solutions whose running maximum of Mercury eccentricity is smaller than
0.7 [7]. Figure 10 shows that the QIs are indeed slow quantities when compared
to the LL action variables. The growth of the QI dispersion is detailed in
Fig. 11, where we report the time evolution of the interquartile range (IQR)
of their distributions. After a transient phase lasting about 100 Myr and
characterized by the exponential separation of close trajectories, the time
growth of the IQR follows a power law typical of diffusion processes. Figures
10 and 11 clearly show the slower diffusion of
$\mathcal{C}_{\textnormal{\scriptsize inc}}$ and $\mathcal{C}_{2}$ in the
model $\mathbb{H}_{2n}^{\bullet}$ when compared to $\mathcal{H}_{2n}$. We
recall that $\mathcal{E}_{2n}^{\bullet}$ is an exact integral of motion for
the model $\mathbb{H}_{2n}^{\bullet}$ (see Sect. IV.4) and its PDF has null
dispersion.

Figure 10:  Time evolution over 100 Gyr of the PDF of the signed deviation
from the mean of the low-pass filtered dimensionless QIs and dimensionless
actions $\mathfrak{X}_{1},\Psi_{3}$. Estimation from an ensemble of 1080
numerical orbital solutions for different models ($\mathcal{H}_{4}$,
$\mathcal{H}_{6}$, $\mathbb{H}_{4}^{\bullet}$, and
$\mathbb{H}_{6}^{\bullet}$). First row:
$\widehat{\mathcal{C}}_{\textnormal{\scriptsize inc}}$. Second row:
$\widehat{\mathcal{C}}_{2}$. Third row: $\widehat{\mathcal{E}}_{4}$
($\mathcal{H}_{4}$) and $\widehat{\mathcal{E}}_{6}$ ($\mathcal{H}_{6}$).
Fourth row: $\widehat{\mathfrak{X}}_{1}$. Fifth row: $\widehat{\Psi}_{3}$. The
time of each curve is color coded. At each time, the estimation only takes
into account stable solutions, that are those with a running maximum of
Mercury eccentricity smaller than 0.7. The quantity
$\widehat{\mathcal{E}}_{2n}^{\bullet}$ is an exact integral of motion for the
model $\mathbb{H}_{2n}^{\bullet}$ and its PDF has null dispersion. Figure 11:
Time evolution of the interquartile range (IQR) of the ensemble PDFs of the
QIs shown in Fig. 10. Left: $\widehat{\mathcal{C}}_{\textnormal{\scriptsize
inc}}$. Middle: $\widehat{\mathcal{C}}_{2}$. Right:
$\widehat{\mathcal{E}}_{4}$ ($\mathcal{H}_{4}$) and
$\widehat{\mathcal{E}}_{6}$ ($\mathcal{H}_{6}$). The quantity
$\widehat{\mathcal{E}}_{2n}^{\bullet}$ is an exact integral of motion for the
model $\mathbb{H}_{2n}^{\bullet}$ and its PDF has a null IQR.

## References

  * Laskar [1989] J. Laskar, A numerical experiment on the chaotic behaviour of the solar system, Nature (London) 338, 237 (1989).
  * Laskar [1990a] J. Laskar, The chaotic motion of the solar system - A numerical estimate of the size of the chaotic zones, Icarus 88, 266 (1990a).
  * Sussman and Wisdom [1992] G. J. Sussman and J. Wisdom, Chaotic evolution of the solar system, Science 257, 56 (1992).
  * Mogavero and Laskar [2021] F. Mogavero and J. Laskar, Long-term dynamics of the inner planets in the Solar System, Astronomy and Astrophysics 655, A1 (2021), arXiv:2105.14976 [astro-ph.EP] .
  * Laskar [2008] J. Laskar, Chaotic diffusion in the Solar System, Icarus 196, 1 (2008), arXiv:0802.3371 .
  * Laskar and Gastineau [2009] J. Laskar and M. Gastineau, Existence of collisional trajectories of Mercury, Mars and Venus with the Earth, Nature (London) 459, 817 (2009).
  * Hoang  _et al._ [2022] N. H. Hoang, F. Mogavero, and J. Laskar, Long-term instability of the inner Solar system: numerical experiments, Monthly Notices of the RAS 514, 1342 (2022), arXiv:2205.04170 [astro-ph.EP] .
  * Batygin  _et al._ [2015] K. Batygin, A. Morbidelli, and M. J. Holman, Chaotic Disintegration of the Inner Solar System, Astrophys. J. 799, 120 (2015), arXiv:1411.5066 [astro-ph.EP] .
  * Woillez and Bouchet [2020] E. Woillez and F. Bouchet, Instantons for the Destabilization of the Inner Solar System, Phys. Rev. Lett. 125, 021101 (2020), arXiv:1910.04005 [nlin.CD] .
  * Note [1] The dynamics truncated at degree 4 produces nevertheless the same chaos of the full system, as measured by the finite-time maximum Lyapunov exponent [19].
  * Milani and Nobili [1992] A. Milani and A. M. Nobili, An example of stable chaos in the Solar System, Nature (London) 357, 569 (1992).
  * Morbidelli and Froeschlé [1996] A. Morbidelli and C. Froeschlé, On the Relationship Between Lyapunov Times and Macroscopic Instability Times, Celestial Mechanics and Dynamical Astronomy 63, 227 (1996).
  * Fermi  _et al._ [1955] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou,  _Studies of Nonlinear Problems_, Tech. Rep. (Los Alamos National Laboratory, 1955).
  * Lam and Kurchan [2014] K.-D. N. T. Lam and J. Kurchan, Stochastic Perturbation of Integrable Systems: A Window to Weakly Chaotic Systems, Journal of Statistical Physics 156, 619 (2014).
  * Goldfriend and Kurchan [2019] T. Goldfriend and J. Kurchan, Equilibration of quasi-integrable systems, Phys. Rev. E 99, 022146 (2019), arXiv:1810.06121 [cond-mat.stat-mech] .
  * Note [2] This resonance involves the fundamental precession frequencies of Mercury and Jupiter perihelia [5, 36, 62].
  * Laskar [1996] J. Laskar, Large Scale Chaos and Marginal Stability in the Solar System, Celestial Mechanics and Dynamical Astronomy 64, 115 (1996).
  * Note [3] “At each stage of its evolution, the system should have a time of stability comparable with its age” [17].
  * Mogavero and Laskar [2022] F. Mogavero and J. Laskar, The origin of chaos in the Solar System through computer algebra, Astronomy and Astrophysics 662, L3 (2022), arXiv:2205.03298 [astro-ph.EP] .
  * Laskar  _et al._ [2004] J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A. C. M. Correia, and B. Levrard, A long-term numerical solution for the insolation quantities of the Earth, Astronomy and Astrophysics 428, 261 (2004).
  * Morbidelli [2002] A. Morbidelli,  _Modern celestial mechanics: aspects of solar system dynamics_ (Taylor & Francis, 2002).
  * Gauss [1818] C. F. Gauss,  _Determinatio attrationis, quam in punctum quodvis positionis datae exerceret planeta, SI eius massa per tota orbitam, ratione temporis, quo singulae partes describuntur, uniformiter errat dispertita_ (Heinrich Dieterich, Göttingen, 1818).
  * Laskar [1991] J. Laskar, Analytical framework in Poincaré variables for the motion of the solar system, in  _Predictability, Stability, and Chaos in N-Body Dynamical Systems_, NATO Advanced Study Institute (ASI) Series B, Vol. 272, edited by A. Roy (Plenum Press, New York, 1991) pp. 93–114.
  * Laskar and Robutel [1995] J. Laskar and P. Robutel, Stability of the Planetary Three-Body Problem. I. Expansion of the Planetary Hamiltonian, Celestial Mechanics and Dynamical Astronomy 62, 193 (1995).
  * Gastineau and Laskar [2011] M. Gastineau and J. Laskar, Trip: A computer algebra system dedicated to celestial mechanics and perturbation series, ACM Commun. Comput. Algebra 44, 194 (2011).
  * Gastineau and Laskar [2021] M. Gastineau and J. Laskar,  _TRIP 1.4.120_, TRIP Reference manual (IMCCE, Paris Observatory, 2021).
  * Note [4] $\mathrm{E}$ represents the exponential operator and $j$ stands for the imaginary unit. The overline on variables denotes complex conjugate.
  * Hoang  _et al._ [2021] N. H. Hoang, F. Mogavero, and J. Laskar, Chaotic diffusion of the fundamental frequencies in the Solar System, Astronomy and Astrophysics 654, A156 (2021), arXiv:2106.00584 [astro-ph.EP] .
  * Laskar [1988] J. Laskar, Secular evolution of the solar system over 10 million years, Astronomy and Astrophysics 198, 341 (1988).
  * Laskar [2005] J. Laskar, Frequency map analysis and quasiperiodic decompositions, in  _Hamiltonian Systems and Fourier Analysis: New Prospects For Gravitational Dynamics_, edited by D. Benest, C. Froeschlé, and E. Lega (Cambridge Scientific Publishers Ltd, Cambridge, 2005) pp. 93–114, arXiv:math/0305364 .
  * Oseledec [1968] V. I. Oseledec, A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems, Trans.Moscow Math.Soc. 19, 197 (1968).
  * Eckmann and Ruelle [1985] J. P. Eckmann and D. Ruelle, Ergodic theory of chaos and strange attractors, Reviews of Modern Physics 57, 617 (1985).
  * Gaspard [1998] P. Gaspard,  _Chaos, Scattering and Statistical Mechanics_, Cambridge Nonlinear Science Series (Cambridge University Press, 1998).
  * Skokos [2010] C. Skokos, The Lyapunov characteristic exponents and their computation, in  _Dynamics of Small Solar System Bodies and Exoplanets_, edited by J. J. Souchay and R. Dvorak (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010) p. 63–135.
  * Laskar [1994] J. Laskar, Large-scale chaos in the solar system., Astronomy and Astrophysics 287, L9 (1994).
  * Batygin and Laughlin [2008] K. Batygin and G. Laughlin, On the Dynamical Stability of the Solar System, Astrophys. J. 683, 1207 (2008), arXiv:0804.1946 [astro-ph] .
  * Note [5] The FT-LCEs do not depend on any tangent vector $\bm{\zeta}_{0}$.
  * Note [6] All the initial conditions mentioned in this work are sampled from a multidimensional Gaussian distribution that is centered at the nominal initial conditions of Gauss’s dynamics [4, Appendix D] and has a relative width of $10^{-9}$ (see [19, Appendix C] and [7, Sect. 3] for details).
  * Morbidelli and Giorgilli [1995] A. Morbidelli and A. Giorgilli, On a connection between KAM and Nekhoroshev’s theorems, Physica D Nonlinear Phenomena 86, 514 (1995).
  * Laskar [1992] J. Laskar, A Few Points on the Stability of the Solar System (lecture), in  _Chaos, Resonance, and Collective Dynamical Phenomena in the Solar System_ , IAU Symposium, Vol. 152, edited by S. Ferraz-Mello (Kluwer Academic Publishers, Dordrecht, 1992) p. 1.
  * Brown and Rein [2020] G. Brown and H. Rein, A Repository of Vanilla Long-term Integrations of the Solar System, Research Notes of the American Astronomical Society 4, 221 (2020), arXiv:2012.05177 [astro-ph.EP] .
  * Benettin  _et al._ [1978] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, All Lyapunov characteristic numbers are effectively computable., Academie des Sciences Paris Comptes Rendus Serie B Sciences Physiques 286, 431 (1978).
  * Benettin  _et al._ [1980] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems - A method for computing all of them. I - Theory. II - Numerical application, Meccanica 15, 9 (1980).
  * Chirikov [1979] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Physics Reports 52, 263 (1979).
  * Laskar [1990b] J. Laskar, Systèmes de Variables et Eléments, in  _Modern Methods in Celestial Mechanics_ , edited by D. Benest and C. Froeschle (Editions Frontières, Gif -Sur-Yvette, 1990) pp. 63–87.
  * Boué and Laskar [2009] G. Boué and J. Laskar, Spin axis evolution of two interacting bodies, Icarus 201, 750 (2009).
  * Laskar [1997] J. Laskar, Large scale chaos and the spacing of the inner planets., Astronomy and Astrophysics 317, L75 (1997).
  * Zurbenko and Smith [2018] I. G. Zurbenko and D. Smith, Kolmogorov–Zurbenko filters in spatiotemporal analysis, WIREs Computational Statistics 10, e1419 (2018).
  * Correia and Laskar [2003] A. C. M. Correia and J. Laskar, Long-term evolution of the spin of Venus. II. numerical simulations, Icarus 163, 24 (2003).
  * Pearson [1901] K. Pearson, LIII. On lines and planes of closest fit to systems of points in space, The London, Edinburgh, and Dublin philosophical magazine and journal of science 2, 559 (1901).
  * Hotelling [1933] H. Hotelling, Analysis of a complex of statistical variables into principal components., Journal of educational psychology 24, 417 (1933).
  * Jolliffe [2002] I. T. Jolliffe,  _Principal component analysis_ (Springer, 2002).
  * Jolliffe and Cadima [2016] I. T. Jolliffe and J. Cadima, Principal component analysis: a review and recent developments, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150202 (2016).
  * Note [7] The projection of a vector $\bm{q}$ onto the subspace spanned by a set of vectors $S$ can be written in a vectorial form as $\mathbf{B}(\mathbf{B}^{\mathrm{T}}\mathbf{B})^{-1}\mathbf{B}^{\mathrm{T}}\bm{q}$, where the column space of the matrix $\mathbf{B}$ is the subspace spanned by the set $S$.
  * Jordan [1875] C. Jordan, Essai sur la géométrie à $n$ dimensions, Bulletin de la Société mathématique de France 3, 103 (1875).
  * Van Loan and Golub [1996] C. F. Van Loan and G. Golub,  _Matrix computations_ (The Johns Hopkins University Press, London, 1996).
  * Ye and Lim [2016] K. Ye and L.-H. Lim, Schubert varieties and distances between subspaces of different dimensions, SIAM Journal on Matrix Analysis and Applications 37, 1176 (2016).
  * Björck and Golub [1973] A. Björck and G. H. Golub, Numerical methods for computing angles between linear subspaces, Mathematics of computation 27, 579 (1973).
  * Knyazev and Argentati [2002] A. V. Knyazev and M. E. Argentati, Principal angles between subspaces in an a-based scalar product: algorithms and perturbation estimates, SIAM Journal on Scientific Computing 23, 2008 (2002).
  * Muller [1959] M. E. Muller, A note on a method for generating points uniformly on n-dimensional spheres, Commun. ACM 2, 19–20 (1959).
  * Mogavero [2017] F. Mogavero, Addressing the statistical mechanics of planet orbits in the solar system, Astronomy and Astrophysics 606, A79 (2017), arXiv:1703.09225 [astro-ph.EP] .
  * Boué  _et al._ [2012] G. Boué, J. Laskar, and F. Farago, A simple model of the chaotic eccentricity of Mercury, Astronomy and Astrophysics 548, A43 (2012), arXiv:1210.5221 [astro-ph.EP] .
  * Hénon [1974] M. Hénon, Integrals of the Toda lattice, Phys. Rev. B 9, 1921 (1974).
  * Flaschka [1974] H. Flaschka, The Toda lattice. II. Existence of integrals, Phys. Rev. B 9, 1924 (1974).
  * Manakov [1974] S. V. Manakov, Complete integrability and stochastization of discrete dynamical systems, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki 67, 543 (1974).
  * Ferguson  _et al._ [1982] J. Ferguson, W. E., H. Flaschka, and D. W. McLaughlin, Nonlinear Normal Modes for the Toda Chain, Journal of Computational Physics 45, 157 (1982).
  * Benettin  _et al._ [2013] G. Benettin, H. Christodoulidi, and A. Ponno, The Fermi-Pasta-Ulam Problem and Its Underlying Integrable Dynamics, Journal of Statistical Physics 152, 195 (2013).
  * Christodoulidi and Efthymiopoulos [2019] H. Christodoulidi and C. Efthymiopoulos, Stages of dynamics in the Fermi-Pasta-Ulam system as probed by the first Toda integral, Mathematics in Engineering 1, 359 (2019).
  * Grava  _et al._ [2020] T. Grava, A. Maspero, G. Mazzuca, and A. Ponno, Adiabatic Invariants for the FPUT and Toda Chain in the Thermodynamic Limit, Communications in Mathematical Physics 380, 811 (2020), arXiv:2001.08070 [math-ph] .
  * Ford [1961] J. Ford, Equipartition of Energy for Nonlinear Systems, Journal of Mathematical Physics 2, 387 (1961).
  * Onorato  _et al._ [2015] M. Onorato, L. Vozella, D. Proment, and Y. V. Lvov, Route to thermalization in the $\alpha$-Fermi-Pasta-Ulam system, Proceedings of the National Academy of Science 112, 4208 (2015), arXiv:1402.1603 [nlin.CD] .
  * Pistone  _et al._ [2019] L. Pistone, S. Chibbaro, M. D. Bustamante, Y. V. Lvov, and M. Onorato, Universal route to thermalization in weakly-nonlinear one-dimensional chains, Mathematics in Engineering 1, 672–698 (2019).
  * Malishava and Flach [2022] M. Malishava and S. Flach, Lyapunov Spectrum Scaling for Classical Many-Body Dynamics Close to Integrability, Phys. Rev. Lett. 128, 134102 (2022), arXiv:2109.01361 [nlin.CD] .
  * Contopoulos [1960] G. Contopoulos, A third Integral of Motion in a Galaxy, Zeitschrift fuer Astrophysik 49, 273 (1960).
  * Kruskal [1962] M. Kruskal, Asymptotic Theory of Hamiltonian and other Systems with all Solutions Nearly Periodic, Journal of Mathematical Physics 3, 806 (1962).
  * Note [8] In this regard, it should be noted that a relation between QIs and Lyapunov exponents has already been highlighted in simple systems [83, 43].
  * Candès  _et al._ [2011] E. J. Candès, X. Li, Y. Ma, and J. Wright, Robust principal component analysis?, Journal of the ACM (JACM) 58, 1 (2011).
  * Markopoulos  _et al._ [2014] P. P. Markopoulos, G. N. Karystinos, and D. A. Pados, Optimal algorithms for $l\\_${$1$}-subspace signal processing, IEEE Transactions on Signal Processing 62, 5046 (2014).
  * Lee and Verleysen [2007] J. A. Lee and M. Verleysen,  _Nonlinear dimensionality reduction_ , Vol. 1 (Springer, 2007).
  * Brumberg and Chapront [1973] V. A. Brumberg and J. Chapront, Construction of a General Planetary Theory of the First Order, Celestial Mechanics 8, 335 (1973).
  * Bretagnon [1974] P. Bretagnon, Long-period terms in the solar system, Astronomy and Astrophysics 30, 141 (1974).
  * Laskar [1985] J. Laskar, Accurate methods in general planetary theory, Astronomy and Astrophysics 144, 133 (1985).
  * Contopoulos  _et al._ [1978] G. Contopoulos, L. Galgani, and A. Giorgilli, On the number of isolating integrals in Hamiltonian systems, Phys. Rev. A 18, 1183 (1978).


# KAM theorem with large perturbation and application to network of Duffing
oscillators

Xiaoping Yuan<EMAIL_ADDRESS>Lu Chen<EMAIL_ADDRESS>Jing Li
<EMAIL_ADDRESS>
###### Abstract

We prove that there is an invariant torus with given Diophantine frequency
vector for a class of Hamiltonian systems defined by an integrable large
Hamiltonian function with a large non-autonomous Hamiltonian perturbation. As
for application, we prove that a finite network of Duffing oscillators with
periodic exterior forces possesses Lagrangian stability for almost all initial
data.

###### keywords:

KAM theorem; Hamiltonian system; invariant torus; Lagrangian stability [2010]
Primary 34L15; Secondary 34B10, 47E05.

††journal: arXiv

[cor]Corresponding author.

##  1 Introduction

It is well-known that the forced harmonic oscillator

$\frac{d^{2}x}{dt^{2}}+\alpha x=k\sin\omega t,\;\alpha>0$ (1.1)

has a periodic solution

$\frac{k}{\alpha-\omega^{2}}\sin\omega t,\;\;(\alpha\neq\omega^{2}).$

When $\alpha=\omega^{2},$ resonance happens. Duffing Duffing  introduced a
nonlinear oscillator with a cubic stiffness term to describe the hardening
spring effect observed in many mechanical problems

$\frac{d^{2}x}{dt^{2}}+\alpha x+\beta x^{2}+\gamma x^{3}=k\sin\omega
t,\;\;\alpha>0,\;\gamma<0.$ (1.2)

Moon-Holmes Moon-Holmes1979 ; Moon-Holmes1980  showed that the equation of the
form

$\frac{d^{2}x}{dt^{2}}+\delta\frac{dx}{dt}-x+x^{3}=k\cos\omega t$ (1.3)

provides the simplest possible model for the forced vibrations of a cantilever
beam in the nonuniform field of two permanent magnets. Such kind of equations
are now called Duffing equations (oscillators). After proving his famous twist
theorem as a part of KAM theory, Moser Moser1962  proposed the study of the
Lagrangian stability for Duffing equation

$\frac{d^{2}x}{dt^{2}}+\alpha x+\beta
x^{3}=p(t),\;\;\alpha\in\mathbb{R},\;\beta>0,$ (1.4)

where $p(t)$ is a periodic function. The Lagrangian stability means that all
solutions exist globally and are bounded in some norm for all time. The proof
of the stability for (1.4) with $\alpha=0,\beta=2$ was completed by Morris
Morris . Subsequently, the Lagrangian stability result was extended by
Dieckerhoff-Zehnder Dieckerhoff1987  to a wider class of systems

$\frac{d^{2}x}{dt^{2}}+x^{2n+1}+\sum_{j=0}^{2n}p_{j}(t)x^{j}=0,\;\;n\geq 1,$
(1.5)

where $p_{j}(t)$’s are smooth enough and of period, say, $2\pi$. See LeviCMP ;
Liu1989 ; Liu1992 ; Yuan1995 ; Yuan1998  for more details.

In many research fields such as physics, mechanics and mathematical biology as
so on arise networks of coupled Duffing oscillators of various form. For
example, the evolution equations for the voltage variables $V_{1}$ and $V_{2}$
obtained using the Kirchhoff’s voltage law are

$\left\\{\begin{array}[]{ll}R^{2}C^{2}\frac{d^{2}V_{1}}{dt^{2}}=-(\frac{R^{2}C}{R_{1}})\frac{dV_{1}}{dt}-(\frac{R}{R_{2}})V_{1}-(\frac{R}{100R_{3}})V_{1}^{3}+(\frac{R}{R_{C}})V_{2}+f\sin\omega
t,\\\
R^{2}C^{2}\frac{d^{2}V_{2}}{dt^{2}}=-(\frac{R^{2}C}{R_{1}})\frac{dV_{2}}{dt}-(\frac{R}{R_{2}})V_{2}-(\frac{R}{100R_{3}})V_{2}^{3}+(\frac{R}{R_{C}})V_{1}\end{array}\right.$
(1.6)

where $R$’s and $C$’s are resistors and capacitors, respectively. This
equation can be regarded as one coupled by two Duffing oscillators. See
Jothimurugan2016 ; Kiss ; Kovaleva ; Clerc2018 ; Sarkar ; Chatterjee ; Shena ;
Sarma  for more details.

In the present paper, we concern the Lagrangian stability for coupled
Hamiltonian system of $m$ Duffing oscillators:

$\ddot{x_{i}}+x_{i}^{2n+1}+\frac{\partial F}{\partial
x_{i}}=0,\;\;i=1,2,\cdots,m,$ (1.7)

where the polynomial potential
$F=F(x,t)=\sum_{\alpha\in\mathbb{N}^{m},|\alpha|\leq
2n+1}p_{\alpha}(t)x^{\alpha},$ $x\in\mathbb{R}^{m}$ with $p_{\alpha}(t)$ is of
period $2\pi$, and $n$ is a given natural number. When $m=1$, (1.7) is reduced
to (1.5). Note that (1.7) is a Hamiltonian system with freedom $m+1/2$ where
time takes $1/2$-freedom. In generally speaking, the study of the stability of
(1.7) with $m>1$ is a difficult task. In order to see this point, let us
recall the method used by Dieckerhoff-Zehnder Dieckerhoff1987  for the case
$m=1$. The key ingredient in Dieckerhoff1987  is the Moser’s twist theorem
Moser1962  which shows that the twist mapping

$P_{\varepsilon}:\left\\{\begin{array}[]{ll}\rho_{1}=\rho+\varepsilon
g(\rho,\theta),&\rho\in[a,b],\\\ \theta_{1}=\theta+\alpha(\rho)+\varepsilon
f(\rho,\theta),&\theta\in\mathbb{T}=\mathbb{R}/2\pi\end{array}\right.$ (1.8)

possesses an invariant curve with rotational frequency
$\omega=\alpha(\rho_{0})\in$ D.C. (D. C. means Diophantine conditions)
provided that $\parallel f\parallel_{C^{333}},$ $\parallel
g\parallel_{C^{333}}\lesssim 1$, $0<\varepsilon\ll 1$ and $P_{\varepsilon}$ is
symplectic or preserving-intersection. ($C^{333}$ can be relaxed to
$C^{3+\mu}$ with $\mu>0$. See Rüssman Russmann .) Observe that the mapping
$P_{\varepsilon}$ is an integrable one attached by a small perturbation
$\varepsilon(f,g)$. Although $\ddot{x}+x^{2n+1}=0$ is indeed integrable, the
perturbation $\sum_{j=0}^{2n}p_{j}(t)x^{j}$ is not small at all in (1.5). An
important observation by Dieckerhoff-Zehnder Dieckerhoff1987  is that the
perturbation $\sum_{j=0}^{2n}p_{j}(t)x^{j}$ is relatively small with respect
to the integrable part $x^{2n+1}$ in the neighborhood of the infinity. More
exactly, one can write (1.5) as a Hamiltonian equation with Hamiltonian $H$:

$H=H_{0}(I)+R(I,\theta,t),\;\;(I,\theta,t)\in\mathbb{R}\times\mathbb{T}\times\mathbb{T},$

where $(I,\theta)$ are action-angle variables, and

$H_{0}(I)\sim I^{\frac{2n+2}{n+2}},\;R(I,\theta,t)\sim I^{\frac{2n+1}{n+2}}.$
(1.9)

Note

$\lim_{I\rightarrow\infty}(H_{0}(I))^{-1}R(I,\theta,t)=0.$ (1.10)

If the integrable Hamiltonian $H_{0}$ and the perturbation $R$ obey (1.10), we
call that $R$ is relatively small with respect to $H_{0}$ in the neighborhood
of the infinity. Then by a series of symplectic transformations, the
relatively small $R$ can be changed into a truly small perturbation. The
symplectic transformations can be implicitly defined by

$\Psi:\left\\{\begin{array}[]{ll}I=\mu+\frac{\partial S}{\partial\theta},\\\
\phi=\theta+\frac{\partial S}{\partial\mu},\end{array}\right.$ (1.11)

where $S=S(\mu,\theta,t)$ is generating function which obeys the homological
equation

$H^{\prime}_{0}(\mu)\frac{\partial
S}{\partial\theta}+R(\mu,\theta,t)=[R](\mu,t),$ (1.12)

where $[R](\mu,t)=\frac{1}{2\pi}\int_{0}^{2\pi}R(\mu,\theta,t)d\theta.$ When
$m=1$, the homological equation is a scalar equation. Thus,

$S=\int_{0}^{\theta}\frac{[R](\mu,t)-R(\mu,\theta,t)}{H^{\prime}_{0}(\mu)}d\theta.$
(1.13)

By (1.9), we have

$S\sim I^{\frac{n+1}{n+2}}.$ (1.14)

Now the perturbation $R$ is changed into $R_{1}$:

$R_{1}=\frac{\partial S}{\partial
t}+\int_{0}^{1}(1-\tau)H^{\prime\prime}_{0}(\mu+\tau\frac{\partial
S}{\partial\theta})(\frac{\partial S}{\partial\theta})^{2}d\tau\sim
O(I^{\frac{n+1}{n+2}})+O(I^{\frac{2n}{n+2}})=O(I^{\frac{2n}{n+2}}).$ (1.15)

Repeat the procedure as the above $\nu$ times with
$\nu>1+\frac{4}{n}+\log_{2}{n}$. Then $R$ is changed into

$R_{\nu}=O(I^{-\frac{c}{n+2}}),\;\;c>0.$ (1.16)

Now $R_{\nu}$ is small when $|I|\rightarrow+\infty.$ By Moser’s twist theorem,
one can prove that there are many invariant cylinder of time period $2\pi$
around $\infty$ in the extended phase $(x,\dot{x},t).$ The solution
$(x(t),\dot{x}(t))$ with the initial data $(x(0),\dot{x}(0))$ in the cylinder
is confined in the cylinder. Thus, the Duffing oscillator (1.5) has Lagrange
stability

$\sup_{t\in\mathbb{R}}(|x(t)|+|\dot{x}(t)|)<\infty.$ (1.17)

When $m\geq 2,$ equation (1.13) does not hold any more. In order to see this
clearly, passing (1.12) to Fourier coefficients of $S$ and $R,$ we have

$\widehat{S}(\mu,k,t)=\frac{{\mathbf{i}}}{\langle\partial_{\mu}\,H_{0}(\mu),k\rangle}\widehat{R}(\mu,k,t),\;\;k\in\mathbb{Z}^{m}\setminus\\{0\\},\,{\mathbf{i}}^{2}=-1,$
(1.18)

where $\widehat{S}(\mu,k,t)$ and $\widehat{R}(\mu,k,t)$ are Fourier
coefficients of $S(\mu,\cdot,t)$ and $R(\mu,\cdot,t),$ respectively. Note that
$0$ is in the closure of the set
$\\{\langle\partial_{\mu}\,H_{0}(\mu),k\rangle:\;k\in\mathbb{Z}^{m}\setminus\\{0\\}\\}$,
that is, $\langle\partial_{\mu}\,H_{0}(\mu),k\rangle$’s are the notorious
small divisors. Thus the estimate (1.14) does not hold when $m\geq 2$.

In the present paper, we directly construct a KAM theorem to deal Hamiltonian
system with large perturbation.

In order to state the theorem, we need some notations. By $C$ (or $c$) denote
a universal positive constant which may be different in different places. When
those constants $C$ and $c$ are necessarily distinguished, denote them by
$C_{0},C_{1},c_{0},c_{1}$, etc. Let positive integer $d$ be the freedom of the
to-be considered Hamiltonian. Let $[1,2]^{d}$ be the product of $d$ intervals
$[1,2]$ and $\mathbb{T}^{d+1}=\mathbb{R}^{d+1}/(2\pi\mathbb{Z})^{d+1}$. Denote
by $\varepsilon>0$ a small constant which measures the size of Hamiltonian
functions.

######  Theorem 1.1.

Consider a Hamiltonian
$H=\varepsilon^{-a}H_{0}(I)+\varepsilon^{-b}R(\theta,t,I)$ where $a,b$ are
given positive constants with $a>b$, and $H_{0}$ and $R$ obey the following
conditions:

  1. 1.

$H_{0}$ and $R(\theta,t,I)$ are real analytic functions in the domain
$\mathbb{T}^{d+1}\times[1,2]^{d}$ and

$||H_{0}||:=\sup_{I\in[1,2]^{d}}|H_{0}(I)|\leq
c_{1},\;||R||:=\sup_{(\theta,t,I)\in\mathbb{T}^{d+1}\times[1,2]^{d}}|R(\theta,t,I)|\leq
c_{2}$

  2. 2.

$H_{0}$ is non-degenerate in Kolmogorov’s sense:

$\text{det}\,\left(\partial_{I}^{2}\,H_{0}(I)\right)\geq
c_{3}>0,\forall\;I\in[1,2]^{d}.$

Then there exists $\epsilon^{*}=\epsilon^{*}(a,b,d,c_{1},c_{2},c_{3})\ll 1$
such that for any $\varepsilon$ with $|\varepsilon|<\epsilon^{*}$, the
Hamiltonian system

$\dot{\theta}=\frac{\partial H(\theta,t,I)}{\partial
I},\;\dot{I}=-\frac{\partial H(\theta,t,I)}{\partial\theta}$

possesses a $d+1$ dimensional invariant torus of rotational frequency vector
$(\omega(I_{0}),2\pi)$ with $\omega(I):=\frac{\partial H_{0}(I)}{\partial I}$,
for any $I_{0}\in[1,2]^{d}$ and $\omega(I_{0})$ obeying Diophantine
conditions:

$|\varepsilon^{-a}\langle
k,\omega(I_{0})\rangle+l|\geq\frac{\varepsilon^{-a}\gamma}{|k|^{d+1}}>\frac{\gamma}{|k|^{d+1}},\;\gamma=(\log\frac{1}{\varepsilon})^{-2C_{1}},k\in\mathbb{Z}^{d}\setminus\\{0\\},l\in\mathbb{Z},|k|+|l|\leq(\log\frac{1}{\varepsilon})^{C_{1}},$

and

$|\varepsilon^{-a}\langle
k,\omega(I_{0})\rangle+l|\geq\frac{\gamma}{(1+|k|)^{d+1}},\;(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z},|k|+|l|>(\log\frac{1}{\varepsilon})^{C_{1}}.$

Applying Theorem 1.1 to (1.7) we have the following theorem.

######  Theorem 1.2.

The equation (1.7) is almost Lagrangian stable, that is, for almost all
initial data $(x_{j}(0),\dot{x}_{j}(0):j=1,...,m)$, the solution
$(x_{j}(t),\dot{x}_{j}(t):j=1,...,m)$ exists globally and obeys

$\sup_{t\in\mathbb{R}}\sum_{i=1}^{m}|x_{i}(t)|+|\dot{x}_{i}(t)|<C$ (1.19)

where $C$ depends on the initial data $(x_{j}(0),\dot{x}_{j}(0):j=1,...,m)$.

######  Remark 1.

Actually we can prove that the solutions to (1.7) are time quasi-periodic for
almost all large initial data. See $\S$ 5.

######  Remark 2.

Theorems 1.1 and 1.2 hold true for the reversible systems, respectively. We do
not pursue this end here.

######  Remark 3.

Let $\varTheta=\\{I_{0}\in[1,2]^{d}:\;\omega(I_{0})\,\text{obeys the
Diophantine conditions}\\}$. We claim that the Lebesgue measure of $\varTheta$
approaches to $1$:

$\text{Leb}\varTheta\geq 1-C(\log\frac{1}{\varepsilon})^{-C}\to
1,\;\text{as}\;\varepsilon\to 0.$

Let

$\tilde{\Theta}_{k,l}=\left\\{\xi\in\omega([1,2]^{d}):\;|\varepsilon^{-a}\langle
k,\xi\rangle+l|\leq\frac{\varepsilon^{-a}\gamma}{|k|^{d+1}},\;0\neq
k\in\mathbb{Z}^{d},\,l\in\mathbb{Z},|k|+|l|\leq(\log\frac{1}{\varepsilon})^{C_{1}}\right\\}$

and

$\tilde{\Theta}_{k,l}=\left\\{\xi\in\omega([1,2]^{d}):\;|\varepsilon^{-a}\langle
k,\xi\rangle+l|\leq\frac{\gamma}{(1+|k|)^{d+1}},\;(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z},|k|+|l|>(\log\frac{1}{\varepsilon})^{C_{1}}\right\\}.$

Note that there exists a direction $k(\xi)$ such that

$\partial_{k(\xi)}(\varepsilon^{-a}\langle
k,\xi\rangle+l)=\varepsilon^{-a}|k|,\;k=(k_{1},...,k_{d})\in\mathbb{Z}^{d}\setminus\\{0\\}.$

Thus,

$\text{Leb}\tilde{\Theta}_{k,l}\leq C\gamma/|k|^{d+2},\;k\neq
0,|k|+|l|\leq(\log\frac{1}{\varepsilon})^{C_{1}},$
$\text{Leb}\tilde{\Theta}_{k,l}\leq C\varepsilon^{a}\gamma/|k|^{d+2},\;k\neq
0,|k|+|l|>(\log\frac{1}{\varepsilon})^{C_{1}}.$

Also note that $\tilde{\Theta}_{k,l}=\emptyset$ when $k=0,l\neq 0$ and when
$|l|>1+\varepsilon^{-a}|k|\sup_{I\in[1,2]^{d}}|\omega(I)|$. Therefore,

$\text{Leb}\,\left(\bigcup_{(0,0)\neq(k,l)\in\mathbb{Z}^{d+1}}\tilde{\Theta}_{k,l}\right)\leq
C\gamma(\log\frac{1}{\varepsilon})^{C_{1}}\leq
C(\log\frac{1}{\varepsilon})^{-C_{1}}.$

Let
$\Theta=[1,2]^{d}\setminus\left(\bigcup_{(0,0)\neq(k,l)\in\mathbb{Z}^{d+1}}\omega^{-1}(\tilde{\Theta}_{k,l})\right)$.
By the Kolmogorov’s non-degenerate condition, the map
$\omega:[1,2]^{d}\to\omega([1,2]^{d})$ is a diffeomorphism in both direction.
Then the proof of the claim is completed by letting
$\Theta=[1,2]^{d}\setminus\left(\bigcup_{(k,l)\in\mathbb{Z}^{d+1}\setminus\\{0,0\\}}\omega^{-1}(\tilde{\Theta}_{k,l})\right)$.

##  2 Normal Form

By the compactness of $\mathbb{T}^{d+1}\times[1,2]^{d}$, there is a constant
$s_{0}>0$ such that the Hamiltonian functions $H_{0}(I),R(\theta,t,I)$ are
analytic in the complex neighborhood
$\mathbb{T}^{d+1}_{2s_{0}}\times[1,2]_{2s_{0}}^{d}$ of
$\mathbb{T}^{d+1}\times[1,2]^{d}$ and are real for real arguments where

$\mathbb{T}^{d+1}_{2s_{0}}\times[1,2]_{2s_{0}}^{d}=\\{\phi\in\mathbb{C}^{d+1}/(2\pi\mathbb{Z})^{d+1}:\;|\Im\phi|\leq
2s_{0}\\}\times\\{z\in\mathbb{C}^{d}:\;\text{dist}_{\mathbb{C}}\,(z,\mathbb{[}1,2]^{d})\leq
2s_{0}\\}.$

We will call that a function of complex variables is real analytic if it is
analytic on the complex variables and real for the real arguments. Thus
$H_{0}$ and $R$ are real analytic and we can assume

$||\partial_{I}^{\alpha}H_{0}||_{2s_{0}}=\sup_{I\in[1,2]^{d}_{2s_{0}}}|\partial_{I}^{\alpha}H_{0}(I)|\leq
C,||\partial^{\alpha}_{I}\partial^{\beta}_{\phi}\,R(\phi,I)||_{2s_{0}}=\sup_{(\phi,I)\in\mathbb{T}^{d+1}_{2s_{0}}\times[1,2]_{2s_{0}}^{d}}|\partial^{\alpha}_{I}\partial^{\beta}_{\phi}R(\phi,I)|\leq
C,\,|\alpha|+|\beta|\leq 2.$ (2.1)

Now take $I_{0}\in[1,2]^{d}$ such that
$\omega(I_{0})=\partial_{I}H_{0}(I_{0})$ obeys D.C.:

$|\varepsilon^{-a}\langle
k,\omega(I_{0})\rangle+l|\geq\frac{\varepsilon^{-a}\gamma}{|k|^{d+1}},\;|k|+|l|\leq{(\log\frac{1}{\varepsilon})^{C_{1}}},\,k\in\mathbb{Z}^{d}\setminus\\{0\\},\,l\in\mathbb{Z}.$
(2.2)

Take a constant $C_{2}$ with $C_{2}>5(d+1)C_{1}$ and define a neighborhood of
$I_{0}$

$B_{0}=\\{I=I_{0}+z:\;z\in\mathbb{C}^{d},\,|z|\leq(\log\frac{1}{\varepsilon})^{-C_{2}}\\}.$

By (2.1) and (2.2), we have that for
$|k|+|l|\leq(\log\frac{1}{\varepsilon})^{C_{1}},\,k\in\mathbb{Z}^{d}\setminus\\{0\\},\,l\in\mathbb{Z},\,I\in
B_{0}$, the frequency $\omega(I)$ obeys

$|\varepsilon^{-a}\langle k,\omega(I)\rangle+l|\geq|\varepsilon^{-a}\langle
k,\omega(I_{0})\rangle+l|-\varepsilon^{-a}|\omega(I)-\omega(I_{0})||k|\geq\frac{\varepsilon^{-a}\gamma}{2|k|^{d+1}}.$
(2.3)

Let $\Omega_{0}=\frac{1}{2}\partial^{2}_{I}H_{0}(I_{0})$. By the Kolmogorov’s
non-degeneracy, we have that $\Omega_{0}$ is invertible matrix and

$||\Omega_{0}||\leq C,\;||\Omega^{-1}_{0}||\leq C.$ (2.4)

Introduce a truncation operator $\Gamma=\Gamma_{K}$ depending on $K>0$ as
follows: For any function $f:\mathbb{T}^{d+1}\rightarrow\mathbb{C}$ (or
$\mathbb{C}^{l},l\geq 1$), write
$f(x)=\sum_{k\in\mathbb{Z}^{d+1}}\widehat{f}(k)e^{i\langle k,x\rangle}$ and
define function $(\Gamma_{K}f)(x)$:

$(\Gamma_{K}f)(x)=\sum_{|k|\leq K}\widehat{f}(k)e^{\mathbf{i}\langle
k,x\rangle}.$

Now consider
$H(\theta,t,I)=\varepsilon^{-a}H_{0}(I)+\varepsilon^{-b}R(\theta,t,I),$ where
$(\theta,t)\in\mathbb{{T}}^{d+1}_{2s_{0}}$ and $I\in B_{0}$. Then $H_{0}(I)$
and $R(\theta,t,I)$ are real analytic on $\mathbb{{T}}^{d+1}_{2s_{0}}\times
B_{0}$. Take $x=(\theta,t).$ Decompose
$R(\theta,t,I)=\Gamma_{K}R+(1-\Gamma_{K})R,$ where

$(\Gamma_{K}R)(x,I)=\sum_{|k|+|l|\leq
K}\widehat{R}(k,l,I)e^{\mathbf{i}(\langle
k,\theta\rangle+l\,t)},\;\;(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z}.$

Let $A=200d(a+b)$ and $m_{0}=2+\left[\frac{b+A}{a-b}\right]$ where $[\cdot]$
is the integer part of a positive number. Define sequences

  * •

$s_{j}=2^{-j}s_{0},\;s_{j}^{l}=s_{j}-\frac{l}{10}(s_{j}-s_{j+1}),l=0,1,...,10,j=0,1,2,...,m_{0};$

  * •

$K_{j}=\frac{A}{s_{j}-s_{j}^{1}}\log\frac{1}{\varepsilon},\;j=0,1,2,...,m_{0};$

  * •

$\tau_{0}=(\log\frac{1}{\varepsilon})^{-C_{2}},\,\tau_{j}=2^{-j}\tau_{0},\tau_{j}^{l}=\tau_{j}-\frac{l}{10}(\tau_{j}-\tau_{j+1}),l=0,1,...,10,\;j=0,1,...,m_{0};$

  * •

$B(\tau_{0})=B_{0},B(\tau_{j})=\\{z\in\mathbb{C}^{d}:\,|z-I_{0}|\leq\tau_{j}\\},\;j=0,1,...,m_{0}$.

Let $C_{1}>1$. We have

$K_{j}\leq K_{m_{0}}<(\log\frac{1}{\varepsilon})^{C_{1}},\;j=0,1,2,...,m_{0}.$
(2.5)

For a function $f$ defined in $\mathbb{T}_{s}^{d+1}\times B(\tau)$ or
$B(\tau)$ , define

$||f||_{s,\tau}=\sup_{(\phi,I)\in\mathbb{T}_{s}^{d+1}\times
B(\tau)}|f(\phi,I)|,\;\text{or}\;||f||_{\tau}=\sup_{I\in B(\tau)}|f(I)|.$

By Cauchy’s formula, we have

$\displaystyle\parallel(1-\Gamma_{K_{0}})R\parallel_{s_{0},\tau_{0}}$
$\displaystyle=$
$\displaystyle\parallel\sum_{|k|+|l|>K_{0}}\widehat{R}(k,l,I)e^{\mathbf{i}(\langle
k,\theta\rangle+l\cdot
t)}\parallel_{s_{0},\tau_{0}}\leq(\sum_{|k|+|l|>K_{0}}e^{-(2s_{0}-s_{0})(|k|+|l|)})\parallel
R\parallel_{2s_{0},2s_{0}}$ $\displaystyle\leq$
$\displaystyle\frac{C}{(s_{0})^{d+1}}e^{-s_{0}\,K_{0}}\leq C\varepsilon^{A},$

and

$\parallel\Gamma_{K_{0}}R\parallel_{s_{0},\tau_{0}}\leq C.$

Let

$H^{(0)}=H,\;\;R^{(0)}=\varepsilon^{-b}\Gamma_{K_{0}}R(\theta,t,I),\;\;R^{(0)}_{+}=\varepsilon^{-b}(1-\Gamma_{K_{0}})R(\theta,t,I).$

Then

$H^{(0)}=\varepsilon^{-a}H_{0}(I)+R^{(0)}(\theta,t,I)+R^{(0)}_{+}(\theta,t,I),\;(\theta,t,I)\in\mathbb{T}^{d+1}_{s_{0}}\times
B(\tau_{0})$ (2.6)

where $R^{(0)}$, $R_{+}^{(0)}$ are real analytic in
$\mathbb{T}^{d+1}_{s_{0}}\times B(\tau_{0}),$ and

$\parallel R^{(0)}\parallel_{s_{0},\tau_{0}}\leq
C\varepsilon^{-b},\Gamma_{K_{0}}\,R^{(0)}\equiv R^{(0)},$ (2.7) $\parallel
R_{+}^{(0)}\parallel_{s_{0},\tau_{0}}\leq C\varepsilon^{A-b}.$ (2.8)

In order to simplify notation, the universal constants $C$’s in the following
lemma may depend on $m_{0},C_{1},C_{2}$ , thus, on $a,b,A,C_{1},C_{2}$.

######  Lemma 2.1 (Iterative lemma in finite steps).

Assume that we have a Hamiltonian

$H^{(j)}=\varepsilon^{-a}H_{0}(I)+h^{(j)}(t,I)+R^{(j)}(\theta,t,I)+R^{(j)}_{+}(\theta,t,I),\,(\theta,t,I)\in\mathbb{T}^{d+1}_{s_{j}}\times
B(\tau_{j}),\,0\leq j\leq m_{0}-1,$ (2.9)

where

1(j)

The functions $h^{(j)}(t,I),R^{(j)}(\theta,t,I)$ and $R^{(j)}_{+}(\theta,t,I)$
are real analytic in $T^{d+1}_{s_{j}}\times B(\tau_{j})$;

2(j)

$h^{(0)}\equiv 0,\;||h^{(j)}||_{s_{j},\tau_{j}}\leq C\varepsilon^{-b},0\leq
j\leq m_{0}-1;$

3(j)

$||R^{(j)}||_{s_{j},\tau_{j}}\leq
C\varepsilon^{-b+j(a-b)}(\log\frac{1}{\varepsilon})^{C},\;{\Gamma_{K_{j}}R^{(j)}=R^{(j)}};$

4(j)

$||R^{(j)}_{+}||_{s_{j},\tau_{j}}\leq
C\varepsilon^{-b+A}(\log\frac{1}{\varepsilon})^{C}.$

Then there exists a symplectic coordinate change

$\Psi_{j}:\;T^{d+1}_{s_{j+1}}\times B(\tau_{j+1})\to T^{d+1}_{s_{j}}\times
B(\tau_{j})$ (2.10)

with

$||\Psi_{j}-id||_{s_{j+1},\tau_{j+1}}\leq
C\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C}$ (2.11)

such that

$H^{(j+1)}=H^{(j)}\circ\Psi_{j}=\varepsilon^{-a}H_{0}(I)+h^{(j+1)}(t,I)+R^{(j+1)}(\theta,t,I)+R^{(j+1)}_{+}(\theta,t,I),\,(\theta,t,I)\in\mathbb{T}^{d+1}_{s_{j+1}}\times
B(\tau_{j+1}),$ (2.12)

where the new Hamiltonian functions $h^{(j+1)}(t,I),R^{(j+1)}(\theta,t,I)$ and
$R^{(j+1)}_{+}(\theta,t,I)$ obey the conditions 1(j+1)-4(j+1).

###### Proof.

Assume that the change $\Psi_{j}$ is implicitly defined by

$\Psi_{j}:\left\\{\begin{array}[]{ll}I=\rho+\frac{\partial
S}{\partial\theta},\\\ \phi=\theta+\frac{\partial S}{\partial\rho},\\\
t=t,\end{array}\right.$ (2.13)

where $S=S(\theta,t,\rho)$ is the generating function, which will be proved to
be analytic in a smaller domain $\mathbb{{T}}^{d+1}_{s_{j+1}}\times
B(\tau_{j+1}).$ By a simple computation, we have

$dI\wedge d\theta=d\rho\wedge
d\theta+\sum_{i,j=1}^{d}\frac{\partial^{2}S}{\partial\rho_{i}\partial\theta_{j}}d\rho_{i}\wedge
d\theta_{j}=d\rho\wedge d\phi.$

Thus the coordinates change $\Psi_{j}$ is symplectic if it exists. Moreover,
we get the changed Hamiltonian

$H^{(j+1)}:=H^{(j)}\circ\Psi_{j}=\varepsilon^{-a}H_{0}(\rho+\frac{\partial
S}{\partial\theta})+h^{(j)}(t,\rho+\frac{\partial
S}{\partial\theta})+R^{(j)}(\theta,t,\rho+\frac{\partial
S}{\partial\theta})+R_{+}^{(j)}(\theta,t,\rho+\frac{\partial
S}{\partial\theta})+\frac{\partial S}{\partial t},$ (2.14)

where $\theta=\theta(\phi,t,\rho)$ is implicitly defined by (2.13). We drop
the sup-index $(j)$. We replace $h^{(j)}$ by $h$, for example.

By Taylor formula, we have

$H^{(j+1)}=\varepsilon^{-a}H_{0}(\rho)+h(t,\rho)+\varepsilon^{-a}\langle\omega(\rho),\frac{\partial
S}{\partial\theta}\rangle+\frac{\partial S}{\partial
t}+R(\theta,t,\rho)+R_{*}+R_{**},$ (2.15)

where

$\displaystyle R_{**}$ $\displaystyle=$ $\displaystyle
R_{+}(\theta,t,\rho+\frac{\partial S}{\partial\theta}),$ (2.16) $\displaystyle
R_{*}$ $\displaystyle=$
$\displaystyle\varepsilon^{-a}\int_{0}^{1}(1-\tau)\partial^{2}_{I}H_{0}(\rho+\tau\frac{\partial
S}{\partial\theta})(\frac{\partial S}{\partial\theta})^{2}d\tau$ (2.19)
$\displaystyle+\int_{0}^{1}\partial_{I}R(\theta,t,\rho+\tau\frac{\partial
S}{\partial\theta})\frac{\partial S}{\partial\theta}d\tau$
$\displaystyle+\int_{0}^{1}\partial_{I}\,h(t,\rho+\tau\frac{\partial
S}{\partial\theta})\frac{\partial S}{\partial\theta}\,d\tau.$

Thus, we derive the homological equation:

$\frac{\partial S}{\partial
t}+\varepsilon^{-a}\langle\omega(\rho),\frac{\partial
S}{\partial\theta}\rangle=\widehat{R}(0,t,\rho)-R(\theta,t,\rho),\quad
S=S(\theta,t,\rho),$ (2.20)

where $\widehat{R}(0,t,\rho)$ is $0$-Fourier coefficient of $R(\theta,t,\rho)$
as the function of $\theta$. Recalling 3(j), we have
$\Gamma_{K_{j}}R^{(j)}=R^{(j)}$. Thus we can assume $\Gamma_{K_{j}}S=S$. For
$f\in\\{S,R\\}$, let

$f(\theta,t,\rho)=\sum_{|k|+|l|\leq K_{j},k\neq
0}\widehat{f}(k,l,\rho)e^{\mathbf{i}(\langle
k,\theta\rangle+lt)},\;\mathbf{i}^{2}=-1.$ (2.21)

By passing to Fourier coefficients, we have

$\widehat{S}(k,l,\rho)=\frac{\mathbf{i}}{\varepsilon^{-a}\langle
k,\omega(\rho)\rangle+l}\widehat{R}(k,l,\rho),\;|k|+|l|\leq
K_{j},\;k\in\mathbb{Z}^{d}\setminus\\{0\\},l\in\mathbb{Z}.$ (2.22)

Let $\rho\in B(\tau_{j})\subset B(\tau_{0})=B_{0}$. By (2.3), we have that

$|\widehat{S}(k,l,\rho)|\leq
C\frac{\varepsilon^{a}|k|^{d+1}}{\gamma}|\widehat{R}(k,l,\rho)|.$

Moreover, by 3(j), for $(\theta,t,\rho)\in\mathbb{T}_{s_{j}^{1}}^{d+1}\times
B(\tau_{j}),$ we get

$|S(\theta,t,\rho)|\leq\frac{C\varepsilon^{a}}{\gamma}\parallel
R\parallel_{s_{j},\tau_{j}}\sum_{|k|+|l|\leq
K_{j}}|k|^{d+1}e^{-(s_{j}-s_{j}^{1})(|k|+|l|)}\leq
C\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C},$ (2.23)

where $C=C(m_{0},d)$. Then by the Cauchy’s estimate, we have

${||\frac{\partial^{p+q}}{\partial\theta^{p}\partial\rho^{q}}S(\theta,t,\rho)||_{s_{j}^{2},\tau_{j}^{1}}\leq
C\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C},p+q\leq 2,p\geq
0,q\geq 0}.$ (2.24)

By (2.13) and (2.24) and the implicit function theorem (analytic version), we
get that there are analytic functions $u=u(\phi,t,\rho),v=v(\phi,t,\rho)$
defined on the domain $\mathbb{T}^{d+1}_{s_{j}^{3}}\times
B(\tau_{j}^{3})\supset\mathbb{T}^{d+1}_{s_{j+1}}\times B(\tau_{j+1})$ with

$\frac{\partial
S(\theta,t,\rho)}{\partial\theta}=u(\phi,t,\rho),\;\frac{\partial
S(\theta,t,\rho)}{\partial\rho}=-v(\phi,t,\rho).$ (2.25)

and

$||u||_{s_{j}^{3},\tau_{j}^{3}}\leq\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C},\;||v||_{s_{j}^{3},\tau_{j}^{3}}\leq\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C}$
(2.26)

such that

$\Psi_{j}:\left\\{\begin{array}[]{ll}I=\rho+u(\phi,t,\rho),\\\
\theta=\phi+v(\phi,t,\rho),\\\ t=t.\end{array}\right.$ (2.27)

Note

$\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C}\ll
s_{j}-s_{j+1},\;\varepsilon^{(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C}\ll\tau_{j}-\tau_{j+1}.$

Then $\Psi_{j}(\mathbb{T}^{d+1}_{s_{j+1}}\times
B(\tau_{j+1}))\subset\mathbb{T}^{d+1}_{s_{j}}\times B(\tau_{j})$. This proves
(2.10) and (2.11). Then by 4(j), we have

$\parallel R_{**}(\phi,t,\rho)\parallel_{s_{j}^{4},\tau_{j}^{4}}\leq\parallel
R_{+}\parallel_{s_{j},\tau_{j}}\leq
C\varepsilon^{-b+A}(\log\frac{1}{\varepsilon})^{C}.$ (2.28)

By (2.25), (2.26), 2(j) and 3(j), we have

$||R_{*}(\phi,t,\rho)||_{s_{j}^{4},\tau_{j}^{4}}\leq
C\left(\varepsilon^{-a}||u||_{s_{j}^{3},\tau_{j}^{3}}^{2}+||R||_{s_{j},\tau_{j}}||u||_{s_{j}^{3},\tau_{j}^{3}}\tau_{j}^{-1}+||h||_{s_{j},\tau_{j}}||u||_{s_{j}^{3},\tau_{j}^{3}}\tau_{j}^{-1}\right)\leq
C\varepsilon^{-b+(j+1)(a-b)}(\log\frac{1}{\varepsilon})^{C}.$ (2.29)

Let

$R^{(j+1)}(\phi,t,\rho)=\Gamma_{K_{j+1}}R_{*}$ (2.30)

and

$R^{(j+1)}_{+}(\phi,t,\rho)=(1-\Gamma_{K_{j+1}})R_{*}+R_{**}.$ (2.31)

By (2.28) and (2.29), we see that $R^{(j+1)}$ and $R^{(j+1)}_{+}$ satisfy
3(j+1) and 4(j+1), respectively. Let

$h^{(j+1)}(t,\rho)=h(t,\rho)+\widehat{R}(0,t,\rho).$ (2.32)

Then by 2(j) and 3(j), we have

$||h^{(j+1)}||_{s_{j+1},\tau_{j+1}}\leq||h^{j}||_{s_{j},\tau_{j}}+||R||_{s_{j},\tau_{j}}\leq
C\varepsilon^{-b},$ (2.33)

which fulfills 2(j+1). Finally, it is obvious that the functions
$h^{(j+1)}(t,I),R^{(j+1)}(\theta,t,I)$ and $R^{(j+1)}_{+}(\theta,t,I)$ are
analytic in $T^{d+1}_{s_{j+1}}\times B(\tau_{j+1})$. Since $R$ is real for the
real arguments, we have

$\widehat{R}(-k,-l,\rho)=\overline{\widehat{R}(k,l,\rho)},\;\rho\in\mathbb{R}^{d},(k,l)\neq(0,0)$

where the bar is the complex conjugate. By (2.22), we have

$\widehat{S}(-k,-l,\rho)=\overline{\widehat{S}(k,l,\rho)},\;\rho\in\mathbb{R}^{d},(k,l)\neq(0,0).$

It follows that $S(\theta,t,\rho)$ is real for the real arguments
$(\theta,t,\rho)$. Moreover, the functions
$h^{(j+1)}(t,I),R^{(j+1)}(\theta,t,I)$ and $R^{(j+1)}_{+}(\theta,t,I)$ are
real for real arguments. By (2.15)-(2.19), (2.30)-(2.32), we have that (2.12)
holds. This completes the proof of the lemma.∎

Applying Lemma 2.1 to (2.6) with $j=0$ and $h^{(0)}\equiv 0$, we get $m_{0}$
symplectic changes $\Psi_{j}$ ($j=0,...,m_{0}-1$) such that

$\Psi:=\Psi_{0}\circ\cdots\circ\Psi_{m_{0}-1},\;\Psi\left(\mathbb{T}^{d+1}_{s_{m_{0}}}\times
B(\tau_{m_{0}})\right)\subset\mathbb{T}^{d+1}_{s_{0}}\times B(\tau_{0}),$
(2.34)
$H^{(m_{0})}=\varepsilon^{-a}H_{0}(\rho)+h^{(m_{0})}(t,\rho)+\tilde{R}(\phi,t,\rho),\,\,(\phi,t,\rho)\in\mathbb{T}^{d+1}_{s_{m_{0}}}\times
B(\tau_{m_{0}}),$ (2.35)

with

$\tilde{R}=\tilde{R}(\phi,t,\rho)=R^{(m_{0})}(\phi,t,\rho)+R^{(m_{0})}_{+}(\phi,t,\rho),\;\,\,(\phi,t,\rho)\in\mathbb{T}^{d+1}_{s_{m_{0}}}\times
B(\tau_{m_{0}}).$ (2.36)

where $h^{(m_{0})},R^{(m_{0})}$ and $R^{(m_{0})}_{+}$ obey the conditions
1(j)-4(j) with j$=m_{0}$. Thus,

$||h^{(m_{0})}(t,\rho)||_{s_{m_{0}},\tau_{m_{0}}}\leq
C\varepsilon^{-b},||\tilde{R}(\phi,t,\rho)||_{s_{m_{0}},\tau_{m_{0}}}\leq
C\varepsilon^{A-b}(\log\frac{1}{\varepsilon})^{C}.$ (2.37)

Let $[h^{(m_{0})}](I)=\widehat{h^{(m_{0})}}(0,I)$ be the $0$-Fourier
coefficient of $h^{(m_{0})}(t,I)$ as the function of $t$. In order to
eliminate the dependence of $h^{(m_{0})}(t,I)$ on the time-variable $t$, we
introduce the following transformation

$\tilde{\Psi}:\;\rho=I,\;\phi=\theta+\frac{\partial}{\partial
I}\tilde{S}(t,I),\;\text{here}\;\tilde{S}(t,I):=\int_{0}^{t}\left([h^{(m_{0})}](I)-h^{(m_{0})}(\xi,I)\right)\,d\xi.$

It is symplectic by easy verification $d\,\rho\wedge d\phi=d\,I\wedge
d\,\theta$. Note that the transformation is not small. So $\tilde{\Psi}$ is
not close to the Identity. In order to apply $\tilde{\Psi}$ to $H^{(m_{0})}$,
we introduce a domain

$\mathcal{D}:=\left\\{t=t_{1}+t_{2}\mathbf{i}\in\mathbb{T}_{s_{m_{0}}}:\;|t_{2}|\leq\varepsilon^{2b}\right\\}\times\left\\{I=x+y\,\mathbf{i}\in
B(\tau_{m_{0}}):\;|x-I_{0}|\leq(\log\frac{1}{\varepsilon})^{-C},\,|y|\leq\varepsilon^{2b}\right\\}.$

Note that $h^{(m_{0})}(t,I)$ is real for real arguments. Thus, for
$(t,I)\in\mathcal{D}$, we have

$\left|\Im\frac{\partial}{\partial
I}\tilde{S}(t,I)\right|=\left|\Im\left(\frac{\partial}{\partial
I}\tilde{S}(t_{1}+t_{2}\mathbf{i},x+y\mathbf{i})-\frac{\partial}{\partial
I}\tilde{S}(t_{1},x)\right)\right|\leq||\partial^{2}\tilde{S}(t,I)||_{\mathcal{D}}(|t_{2}|+|y|)<\varepsilon^{\frac{b}{2}}<\frac{1}{2}s_{m_{0}}.$

Therefore,

$\tilde{\Psi}(\mathbb{T}^{d}_{s_{m_{0}}/2}\times\mathcal{D})\subset\mathbb{T}^{d+1}_{s_{m_{0}}}\times
B(\tau_{m_{0}}),$

and

$\tilde{H}:=H^{(m_{0})}\circ\tilde{\Psi}=\varepsilon^{-a}H_{0}(I)+[h^{(m_{0})}](I)+\breve{R}(\theta,t,I),\;(\theta,t,I)\in\mathbb{T}^{d}_{s_{m_{0}}/2}\times\mathcal{D}.$
(2.38)

where

$\breve{R}(\theta,t,I):=\tilde{R}(\theta+\partial_{I}\tilde{S}(t,I),t,I).$
(2.39)

Recall that $\text{det}\left(\partial_{I}^{2}\,H_{0}(I)\right)\geq c_{3}>0$.
By 2(j) with $j=m_{0}$, we have

$\varepsilon^{a}||\partial_{I}^{2}[h^{(m_{0})}](I)||_{\tau_{m_{0}}/2}\leq
C\varepsilon^{a-b}\left(\log\frac{1}{\varepsilon}\right)^{C}\ll c_{3}.$ (2.40)

Solving the equation
$\partial_{I}H_{0}(I)+\varepsilon^{a}\partial_{I}[h^{(m_{0})}](I)=\omega(I_{0})$
by Newton iteration, we get that there exists $I_{*}\in\mathbb{R}^{d}\bigcap
B(\tau_{m_{0}}/2)$ with $|I_{*}-I_{0}|\leq
C\varepsilon^{a-b}(\log\frac{1}{\varepsilon})^{C}$ such that

$\partial_{I}H_{0}(I_{*})+\varepsilon^{a}\partial_{I}[h^{(m_{0})}](I_{*})=\omega(I_{0}),\;\text{here}\;\omega(I_{0})=\partial_{I}H_{0}(I_{0}).$

Let $I=I_{*}+\rho$ with $\rho\in\mathbb{C}^{d}$ and $|\rho|<\varepsilon^{2b}$.
By Taylor formula at $I=I_{*}$ up to the second order,

$\tilde{H}=\varepsilon^{-a}H_{0}(I_{*})+[h^{(m_{0})}](I_{*})+\varepsilon^{-a}\langle\omega(I_{0}),\rho\rangle+\varepsilon^{-a}\langle\Omega\rho,\rho\rangle+R_{low}(\theta,t,\rho)+R_{high}(\theta,t,\rho)$
(2.41)

where

$\Omega=\frac{1}{2}\partial^{2}_{I}\left.\left(H_{0}(I)+\varepsilon^{a}[h^{(m_{0})}](I)\right)\right|_{I=I_{*}}$
(2.42)
$R_{low}(\theta,t,\rho)=\breve{R}(\theta,t,I_{*})+\langle\partial_{I}\breve{R}(\theta,t,I_{*}),\rho\rangle+\frac{1}{2}\langle\partial_{I}^{2}\breve{R}(\theta,t,I_{*})\rho,\rho\rangle$
(2.43)
$R_{high}(\theta,t,\rho)=\varepsilon^{-a}\int_{0}^{1}\frac{(1-x)^{2}}{2}\partial^{3}_{x}\left(H_{0}(I_{*}+x\rho)+\varepsilon^{a}[h^{(m_{0})}](I_{*}+x\rho)+\varepsilon^{a}\breve{R}(\theta,t,I_{*}+x\rho)\right)\,dx.$
(2.44)

We introduce some notations: With abuse of notation a bit, let
$B(r)=\\{z\in\mathbb{C}^{d}:\;|z-I_{*}|\leq r\\}$. If
$f:\mathbb{T}_{s}^{d+1}\times B(r)\rightarrow\mathbb{C}^{l}(l\geq 1)$ is real
analytic, and

$\sup_{x\in\mathbb{T}^{d+1}_{s}\times B(r)}|f(x,I)|\leq
C\varepsilon^{\beta}|I|^{\alpha},\alpha\geq 0,\beta\in\mathbb{R},$

then we write $f=O_{s,r}(\varepsilon^{\beta}|I|^{\alpha}).$ If $f$ is
independent of $I$, furthermore, we write
$f=O_{s}(\varepsilon^{\beta}|I|^{\alpha}).$ Note that
$\varepsilon(\log\frac{1}{\varepsilon})^{C}\ll 1$. Let $E=A-9b-1$. In view of
(2.37), we have

$\partial_{I}^{p}\breve{R}(\theta,t,I_{*})=O_{\varepsilon^{2b}}(\varepsilon^{E}),p=0,1,2,$
(2.45)
$R_{high}(\theta,t,\rho)=O_{\varepsilon^{2b},\varepsilon^{2b}}(\varepsilon^{-(a+6b)}|\rho|^{3}).$
(2.46)

By (2.4), (2.40) and (2.42), we have

$||\Omega||\leq C,\quad||\Omega^{-1}||\leq C.$ (2.47)

##  3 Extension of Kolmogorov’s Theorem

By a bit abuse of notation, reset $s_{0}:=\varepsilon^{2b}$ and
$r_{0}:=\varepsilon^{2b}$.

######  Theorem 3.1.

Consider a Hamiltonian

$H=C+N(I)+R_{low}(\theta,t,I)+R_{high}(\theta,t,I),\;(\theta,t,I)\in\mathbb{T}^{d+1}_{s_{0}}\times
B(r_{0})$

with the symplectic structure $dI\wedge d\theta$ satisfying

  * (1)

$N(I)=\varepsilon^{-a}\langle\omega,I\rangle+\varepsilon^{-a}\langle\Omega
I,I\rangle,$

  * (2)

$\omega=\omega(I_{0})$ is D.C., that is

$|\varepsilon^{-a}\langle
k,\omega\rangle+l|\geq\frac{\gamma}{(1+|k|)^{d+1}},(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z}\setminus\\{(0,0)\\},$

  * (3)

$||\Omega||\leq C,\;||\Omega^{-1}||\leq C,$

  * (4)

$N(I),$ $R_{low}(\theta,t,I)$ and $R_{high}(\theta,t,I)$ are real analytic in
$\mathbb{T}^{d+1}_{s_{0}}\times B(r_{0}),$

  * (5)

$R_{low}=R_{0}(\theta,t)+\langle R_{1}(\theta,t),I\rangle+\langle
R_{2}(\theta,t)I,I\rangle,$ with $R_{j}=O_{s_{0}}(\varepsilon^{E}),$
$j=0,1,2$,

  * (6)

$R_{high}=O_{s_{0},r_{0}}(\varepsilon^{-(a+6b)}|I|^{3}).$

Then there exists a symplectic coordinate change

$\Psi:\;I=\rho+u(\phi,t,\rho),\theta=\phi+v(\phi,t,\rho),t=t$

where
$u=O_{s_{0}/2,r_{0}/2}(\varepsilon^{\frac{E}{2}}),v=O_{s_{0}/2,r_{0}/2}(\varepsilon^{\frac{E}{2}})$
such that $\Psi\left(\mathbb{T}^{d+1}_{s_{0}/2}\times
B(\frac{r_{0}}{2})\right)\subset\mathbb{T}^{d+1}_{s_{0}}\times B(r_{0})$ and
the Hamiltonian $H$ is changed by $\Psi$ into

$H^{\infty}(\phi,t,\rho)=H\circ\Psi(\phi,t,\rho)=C+\varepsilon^{-a}\langle\omega,\rho\rangle+\varepsilon^{-a}\langle\Omega^{\infty}\rho,\rho\rangle+O_{s_{0}/2,r_{0}/2}(\varepsilon^{-a-6b}|\rho|^{3}),(\phi,t,\rho)\in\mathbb{T}^{d+1}_{s_{0}/2}\times
B(\frac{r_{0}}{2}),$

where

$|\Omega^{\infty}-\Omega|\leq C\varepsilon^{\frac{E}{2}}.$

######  Corollary 3.2.

The Hamiltonian system defined by Hamiltonian $H$ has an invariant torus with
rotational frequency $(\varepsilon^{-a}\omega(I_{0}),2\pi)$ in the extended
phase space $\mathbb{T}^{d+1}\times\mathbb{R}^{d}$.

The proof is finished by the following iterative lemma. To this end, we
introduce some iterative constants, iterative parameters and iterative
domains:

  * •

$\epsilon_{0}=\varepsilon^{E},\epsilon_{j}=(\varepsilon^{E})^{(4/3)^{j}},j=1,2,\cdots,\epsilon_{0}>\epsilon_{1}>\cdots>\epsilon_{j}\downarrow
0$;

  * •

$e_{j}=\frac{\sum_{l=1}^{j}l^{-2}}{100\sum_{l=1}^{\infty}l^{-2}}$;

  * •

$s_{j}=s_{0}(1-e_{j})$ (so $s_{j}>s_{0}/2$ for all $j=1,2,\cdots$);

  * •

$r_{j}=r_{0}(1-e_{j})$ (so $r_{j}>r_{0}/2$ for all $j=1,2,\cdots$);

  * •

$s_{j}^{l}=(1-\frac{l}{10})s_{j}+\frac{l}{10}s_{j+1}$ (so
$s_{j+1}<s_{j}^{l}<s_{j},l=1,2,\cdots,10$);

  * •

$r_{j}^{l}=(1-\frac{l}{10})r_{j}+\frac{l}{10}r_{j+1}$ (so
$r_{j+1}<r_{j}^{l}<r_{j},l=1,2,\cdots,10$);

  * •

For any $\alpha>0$, write
$\epsilon_{j}^{\alpha-}=j^{C}\,\left(\log\frac{1}{\varepsilon}\right)^{C}\,\epsilon^{\alpha}_{j}$.

######  Lemma 3.3 (Iterative Lemma).

Assume that we have had $m$ Hamiltonian functions, for $j=0,1,\cdots,m-1$
satisfying

$\displaystyle H^{(j)}=N^{(j)}+R_{low}^{(j)}+R_{high}^{(j)},$ (3.1)
$\displaystyle
N^{(j)}=\varepsilon^{-a}\langle\omega,I\rangle+\varepsilon^{-a}\langle\Omega^{(j)}I,I\rangle,\parallel\Omega^{(j)}-\Omega^{(j-1)}\parallel\leq
C\epsilon_{j-1}\varepsilon^{-b^{\prime}},j\geq
1,\Omega^{(0)}=\Omega,b^{\prime}:=12d(a+b),$ (3.2) $\displaystyle
R_{low}^{(j)}(\theta,t,I),R_{high}^{(j)}(\theta,t,I)\text{ are real analytic
in }\mathbb{T}^{d+1}_{s_{j}}\times B(r_{j}),$ (3.3) $\displaystyle
R^{(j)}_{low}=R^{(j)}_{0}(\theta,t)+\langle
R^{(j)}_{1}(\theta,t),I\rangle+\langle R^{(j)}_{2}(\theta,t)I,I\rangle,$ (3.4)
$\displaystyle R^{(j)}_{p}(\theta,t,I)=O_{s_{j}}(\epsilon_{j}),\;(p=0,1,2),$
(3.5) $\displaystyle
R_{high}^{(j)}(\theta,t,I)=O_{s_{j},r_{j}}(\varepsilon^{-(a+6b)}|I|^{3}).$
(3.6)

Then there exists a symplectic coordinate

$\Psi_{m}:\;I=\nu+\rho+u_{m}(\phi,t,\rho),\theta=\phi+v_{m}(\phi,t,\rho),t=t$
(3.7)

where $\nu$ is a constant depending $\varepsilon$,
$(\phi,t,\rho)\in\mathbb{T}^{d+1}_{s_{m}}\times B(r_{m})$ and

$u_{m}=O_{s_{m},r_{m}}(\varepsilon^{-3b^{\prime}}\epsilon_{m-1}),v_{m}=O_{s_{m},r_{m}}(\varepsilon^{-3b^{\prime}}\epsilon_{m-1}),$
(3.8) $\Psi_{m}(\mathbb{T}^{d+1}_{s_{m}}\times
B(r_{m}))\subset\mathbb{T}^{d}_{s_{m-1}}\times B(r_{m-1})$ (3.9)

such that the changed Hamiltonian $H^{(m)}=H^{(m-1)}\circ\Psi_{m}$ satisfies
(3.1)-(3.6) by replacing $j$ by $m.$

###### Proof.

Assume the coordinate change $\Psi=\Psi_{m}$ can be implicitly defined by

$\Psi:\;I=\nu+\rho+\frac{\partial
S}{\partial\theta},\phi=\theta+\frac{\partial S}{\partial\rho},t=t,$ (3.10)

where function $S=S(\theta,t,\rho)$ and constant vector $\nu$ will be
specified by homological equations later. Then $\Psi$ is symplectic if it is
well-defined. Let

$S=S_{0}(\theta,t)+\langle S_{1}(\theta,t),\rho\rangle+\langle
S_{2}(\theta,t)\rho,\rho\rangle.$ (3.11)

Then

$\displaystyle H^{(m)}$ $\displaystyle=$ $\displaystyle
H^{(m-1)}\circ\Psi_{m}$ $\displaystyle=$ $\displaystyle
N^{(m-1)}(v+\rho+\frac{\partial
S}{\partial\theta})+R_{low}^{(m-1)}(\theta,t,v+\rho+\frac{\partial
S}{\partial\theta})+R_{high}^{(m-1)}(\theta,t,v+\rho+\frac{\partial
S}{\partial\theta})+\frac{\partial S(\theta,t,\rho)}{\partial t}$
$\displaystyle=$
$\displaystyle\varepsilon^{-a}\langle\omega,v+\rho+\partial_{\theta}S\rangle+\varepsilon^{-a}\langle\Omega^{(m-1)}(v+\rho+\partial_{\theta}S),v+\rho+\partial_{\theta}S\rangle+R_{0}^{(m-1)}(\theta,t)+\frac{\partial
S}{\partial t}$ $\displaystyle+\langle
R_{1}^{(m-1)}(\theta,t),v+\rho+\partial_{\theta}S\rangle+\langle
R_{2}^{(m-1)}(\theta,t)(v+\rho+\partial_{\theta}S),v+\rho+\partial_{\theta}S\rangle+R_{high}(\theta,t,v+\rho+\partial_{\theta}S).$

Let $\omega\cdot\partial_{\theta}S=\langle\omega,\partial_{\theta}S\rangle.$
By a simple computation, we have

$\displaystyle H^{(m)}$ $\displaystyle=$
$\displaystyle\varepsilon^{-a}\langle\omega,v\rangle+\varepsilon^{-a}\langle\omega,\rho\rangle+\varepsilon^{-a}\left(\omega\cdot\partial_{\theta}S_{0}+\langle\omega\cdot\partial_{\theta}S_{1},\rho\rangle+\langle\omega\cdot\partial_{\theta}S_{2}\rho,\rho\rangle\right)$
$\displaystyle+\varepsilon^{-a}\left(\langle\Omega^{(m-1)}v,v\rangle+2\langle\Omega^{(m-1)}v,\rho\rangle+\langle\Omega^{(m-1)}\rho,\rho\rangle+2\langle\Omega^{(m-1)}v,\partial_{\theta}S\rangle\right)$
$\displaystyle+\varepsilon^{-a}\left(2\langle\Omega^{(m-1)}\rho,\partial_{\theta}S_{0}\rangle+2\langle\Omega^{(m-1)}\rho,\langle\partial_{\theta}S_{1},\rho\rangle\rangle+2\langle\Omega^{(m-1)}\rho,\langle\partial_{\theta}S_{2}\rho,\rho\rangle\rangle\right)$
$\displaystyle+\varepsilon^{-a}\langle\Omega^{(m-1)}\partial_{\theta}S,\partial_{\theta}S\rangle+R_{0}^{(m-1)}+\langle
R_{1}^{(m-1)},v+\partial_{\theta}S\rangle+\langle R_{1}^{(m-1)},\rho\rangle$
$\displaystyle+\langle
R_{2}^{(m-1)}(v+\partial_{\theta}S),v+\partial_{\theta}S\rangle+\langle
R_{2}^{(m-1)}\rho,\rho\rangle+2\langle
R_{2}^{(m-1)}(v+\partial_{\theta}S),\rho\rangle$
$\displaystyle+\frac{1}{2}\langle(\partial_{I}^{3}R_{high}(\theta,t,0)\partial_{\theta}S_{0})\rho,\rho\rangle+\partial_{t}S$
$\displaystyle+R_{high}(\theta,t,v+\rho+\partial_{\theta}S)-\frac{1}{2}\langle(\partial_{I}^{3}R_{high}(\theta,t,0)\partial_{\theta}S_{0})\rho,\rho\rangle.$

In the following, we omit the term
$\varepsilon^{-a}(\langle\omega,v\rangle+\langle\Omega^{(m-1)}v,v\rangle)$
which does not affect the dynamics. Let

$R_{*}(\theta,t)=2\varepsilon^{-a}\Omega^{(m-1)}\partial_{\theta}S_{0}(\theta,t)+R_{1}^{(m-1)}(\theta,t),$
(3.13)
$R_{**}(\theta,t)=\frac{1}{2}\partial^{3}_{I}R_{high}(\theta,t,0)\partial_{\theta}S_{0}(\theta,t)+2\varepsilon^{-a}\Omega^{(m-1)}(\partial_{\theta}S_{1}(\theta,t))+R_{2}^{(m-1)}(\theta,t).$
(3.14)

From the changed $H^{(m)}$, we derive the homological equations:

$\displaystyle\varepsilon^{-a}\omega\cdot\partial_{\theta}S_{0}(\theta,t)+\partial_{t}S_{0}(\theta,t)+R_{0}^{(m-1)}(\theta,t)=\widehat{R_{0}^{(m-1)}}(0,0),$
(3.15)
$\displaystyle\varepsilon^{-a}\omega\cdot\partial_{\theta}S_{1}(\theta,t)+\partial_{t}S_{1}(\theta,t)+R_{*}(\theta,t)=\widehat{R_{*}}(0,0),$
(3.16) $\displaystyle 2\varepsilon^{-a}\Omega^{(m-1)}v=-\widehat{R_{*}}(0,0),$
(3.17)
$\displaystyle\varepsilon^{-a}\omega\cdot\partial_{\theta}S_{2}(\theta,t)+\partial_{t}S_{2}(\theta,t)+R_{**}(\theta,t)=\widehat{R_{**}}(0,0).$
(3.18)

Let

$\Omega^{(m)}=\Omega^{(m-1)}+\varepsilon^{a}\widehat{R_{**}}(0,0),$ (3.19)
$\displaystyle R_{+}$ $\displaystyle=$
$\displaystyle\varepsilon^{-a}\left(\langle\Omega^{(m-1)}\partial_{\theta}S,\partial_{\theta}S\rangle+2\langle\Omega^{(m-1)}\nu,\partial_{\theta}S\rangle\right)$
(3.24) $\displaystyle+\langle R_{1}^{(m-1)},\nu+\partial_{\theta}S\rangle$
$\displaystyle+\langle
R_{2}^{(m-1)}(\nu+\partial_{\theta}S),\nu+\partial_{\theta}S\rangle+2\langle
R_{2}^{(m-1)}(\nu+\partial_{\theta}S),\rho\rangle$
$\displaystyle+R_{high}(\theta,t,v+\rho+\partial_{\theta}S)-\frac{1}{2}\langle(\partial_{I}^{3}R_{high}(\theta,t,0)\partial_{\theta}S_{0})\rho,\rho\rangle$
$\displaystyle+2\varepsilon^{-a}\langle\Omega^{(m-1)}\rho,\langle\partial_{\theta}S_{2}\rho,\rho\rangle\rangle.$

Then

$H^{(m)}=C+\varepsilon^{-a}(\langle\omega,\rho\rangle+\varepsilon^{-a}\langle\Omega^{(m)}\rho,\rho\rangle)+R_{+}.$
(3.25)

$\bullet$ Soluitons of the homological equations

  * (1)

Solution to (3.15): Passing to Fourier coefficients,

${S}_{0}(\theta,t)=\sum_{(0,0)\neq(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z}}\frac{\mathbf{i}}{\varepsilon^{-a}\langle
k,\omega\rangle+l}\widehat{R^{(m-1)}_{0}}(k,l)e^{\mathbf{i}(\langle
k,\theta\rangle+lt)}.$

By $\omega\in D.C.,$

$\parallel
S_{0}\parallel_{s_{m-1}^{1}}\leq\frac{C}{(s_{m-1}-s_{m-1}^{1})^{2d+2}}\frac{1}{\gamma}\parallel
R_{0}^{(m-1)}\parallel_{s_{m-1}}\leq\epsilon^{1-}_{m-1}\varepsilon^{-b(4d+4)}.$
(3.26)

  * (2)

Solution to (3.16): By (3.5), (3.13) and (3.26) with $j=m-1$, and using
Cauchy’s estimate, we have

$R_{*}=O_{s_{m-1}^{2}}(\epsilon_{m-1}^{1-}\varepsilon^{-b(4d+6)-a}).$ (3.27)

By (3.16), we have

${S}_{1}(\theta,t)=\sum_{(0,0)\neq(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z}}\frac{\mathbf{i}}{\varepsilon^{-a}\langle
k,\omega\rangle+l}\widehat{R_{*}}(k,l)e^{\mathbf{i}(\langle
k,\theta\rangle+lt)}.$

Therefore

$\parallel
S_{1}\parallel_{s_{m-1}^{3}}\leq\frac{C}{(s_{m-1}^{2}-s_{m-1}^{3})^{2d+2}}\frac{1}{\gamma}\parallel
R_{*}\parallel_{s_{m-1}^{2}}\leq\epsilon^{1-}_{m-1}\varepsilon^{-b(8d+10)-a}.$
(3.28)

  * (3)

Solution to (3.17): By (3.17) and (3.27),

$|\nu|=\left|(2\varepsilon^{-a}\Omega^{(m-1)})^{-1}\widehat{R_{*}}(0,0)\right|\leq
C\epsilon_{m-1}^{1-}\varepsilon^{-b(4d+6)}.$ (3.29)

  * (4)

Similarly, by (3.14), (3.26), (3.28) and Cauchy’s estimate, we have

$\parallel R_{**}(\theta,t)\parallel_{s^{4}_{m-1}}\leq
C\epsilon^{1-}_{m-1}\varepsilon^{-b(8d+12)-2a}.$

Passing (3.18) to Fourier coefficients,

${S}_{2}(\theta,t)=\sum_{(0,0)\neq(k,l)\in\mathbb{Z}^{d}\times\mathbb{Z}}\frac{\mathbf{i}}{\varepsilon^{-a}\langle
k,\omega\rangle+l}\widehat{R_{**}}(k,l)e^{\mathbf{i}(\langle
k,\theta\rangle+lt)}.$

It follows

$\parallel
S_{2}\parallel_{s_{m-1}^{5}}\leq\frac{C}{(s_{m-1}^{4}-s_{m-1}^{5})^{2d+2}}\frac{1}{\gamma}\parallel
R_{**}\parallel_{s_{m-1}^{4}}\leq\epsilon^{1-}_{m-1}\varepsilon^{-b(12d+16)-2a}.$
(3.30)

Now

$|\Omega^{(m)}-\Omega^{(m-1)}|=\varepsilon^{a}|\widehat{R_{**}}(0,0)|\leq\varepsilon^{-b^{\prime}}\epsilon_{m-1}.$
(3.31)

This verifies (3.2) in the Iterative Lemma. By (3.26), (3.28) and (3.30), we
have

$\parallel S\parallel_{s_{m-1}^{5},r_{m-1}}=\parallel S_{0}+\langle
S_{1},\rho\rangle+\langle
S_{2}\rho,\rho\rangle\parallel_{s_{m-1}^{5},r_{m-1}}\leq
C\varepsilon^{-\frac{3b^{\prime}}{2}}\epsilon_{m-1}.$ (3.32)

By Cauchy’s estimate,

$\parallel\sum_{|\alpha|+|\beta|\leq
2}\partial_{\theta}^{\alpha}\partial_{\rho}^{\beta}S\parallel_{s_{m-1}^{6},r_{m-1}^{1}}\leq
C\varepsilon^{-2b^{\prime}}\epsilon_{m-1}.$ (3.33)

By using the implicit theorem, we have that there exist
$u=u_{m}(\phi,t,\rho),$ $v=v_{m}(\phi,t,\rho)$ with

$\parallel u\parallel_{s_{m-1}^{7},r_{m-1}^{7}},\parallel
v\parallel_{s_{m-1}^{7},r_{m-1}^{7}}\leq
C\varepsilon^{-3b^{\prime}}\epsilon_{m-1}$ (3.34)

such that $\Psi=\Psi_{m}$

$\Psi_{m}:\left\\{\begin{array}[]{ll}I=v+\rho+u(\phi,t,\rho),\\\
\theta=\phi+v(\phi,t,\rho),\\\
t=t\end{array}\right.(\phi,t,\rho)\in\mathbb{T}_{s_{m-1}^{7}}^{d+1}\times
B(r_{m-1}^{7})$ (3.35)

is well-defined and

$\Psi_{m}(\mathbb{T}_{s_{m}}^{d+1}\times
B(r_{m}))\subset\mathbb{T}_{s_{m-1}}^{d+1}\times B(r_{m-1}).$

This verifies (3.8) and (3.9).

We are now in position to estimate the new perturbation $R_{+}.$ By (3.29),
(3.33), (3.34) and (3.35), we have

$\displaystyle\parallel\eqref{eq82}\parallel_{s_{m-1}^{7},r_{m-1}^{7}}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-a}(\varepsilon^{-3b^{\prime}}\epsilon_{m-1}^{1-})^{2}\leq
C\varepsilon^{-a-6b^{\prime}}\epsilon_{m-1}^{2}\left(\log\frac{1}{\varepsilon}\right)^{C}\leq\epsilon_{m}\varepsilon^{4b^{\prime}}.$
(3.36)

Similarly,

$\parallel\eqref{eq83}\parallel_{s_{m-1}^{7},r_{m-1}^{7}}\leq\epsilon_{m}\varepsilon^{4b^{\prime}},$
(3.37)
$\parallel\eqref{eq84}\parallel_{s_{m-1}^{7},r_{m-1}^{7}}\leq\epsilon_{m}\varepsilon^{4b^{\prime}}.$
(3.38)

Clearly,

$\eqref{eq86}=O_{s_{m-1}^{7},r_{m-1}^{7}}(\varepsilon^{-4b^{\prime}}\epsilon_{m-1}^{1-}|\rho|^{3})=O_{s_{m-1}^{7},r_{m-1}^{7}}(|\rho|^{3}).$
(3.39)

Recall $R_{high}=O_{s_{m-1},r_{m-1}}(\varepsilon^{-(a+6b)}|I|^{3}).$ Then by
Taylor formula and in view of (3.32), we have

$\displaystyle\eqref{eq85}$ $\displaystyle=$ $\displaystyle
R_{high}(\theta,t,\rho+\partial_{\theta}S)-\langle(\partial_{I}^{3}R_{high}(\theta,t,0)\partial_{\theta}S_{0})\rho,\rho\rangle$
(3.40) $\displaystyle=$ $\displaystyle\sum_{|\alpha|\geq
3,\alpha\in\mathbb{Z}^{d}_{+}}\frac{1}{\alpha!}\,\left((\rho+\partial_{\theta}\,S)\cdot\partial_{I}\right)^{\alpha}R_{high}(\theta,t,0)-\frac{1}{2}\langle(\partial_{I}^{3}R_{high}(\theta,t,0)\partial_{\theta}S_{0})\rho,\rho\rangle$
$\displaystyle=$ $\displaystyle
O_{s_{m-1}^{7},r_{m-1}^{7}}(\varepsilon^{-(a+6b)}|\rho|^{3})+O_{s_{m-1}^{7},r_{m-1}^{7}}(\epsilon_{m}\varepsilon^{3b^{\prime}}).$

Consequently, by (3.36)-(3.40), we have

$R_{+}(\phi,t,\rho)=O_{s_{m-1}^{7},r_{m-1}^{7}}(\varepsilon^{-(a+6b)}|\rho|^{3})+O_{s_{m-1}^{7},r_{m-1}^{7}}(\epsilon_{m}\varepsilon^{3b^{\prime}}).$
(3.41)

By developing $R_{+}(\phi,t,\rho)$ into Taylor formula of order $3$ in $\rho$,
we can re-write

$R_{+}(\phi,t,\rho)=R^{(m)}_{low}(\phi,t,\rho)+R^{(m)}_{high}(\phi,t,\rho)$

where

$R^{(m)}_{low}=R_{0}^{(m)}(\phi,t)+\langle
R_{1}^{(m)}(\phi,t),\rho\rangle+\langle
R_{2}^{(m)}(\phi,t)\,\rho,\rho\rangle,$

and $R_{p}^{(m)}(\phi,t)=O_{s_{m},r_{m}}(\epsilon_{m})$ ($p=0,1,2$) and
$R^{(m)}_{high}(\phi,t,\rho)=O_{s_{m},r_{m}}(\varepsilon^{-(a+6b)}|\rho|^{3})$.
This verifies (3.4)-(3.6) with $j=m$. As for the functions
$R^{(m)}_{low}(\phi,t,\rho)$ and $R^{(m)}_{high}(\phi,t,\rho)$ is real for the
real argument $(\phi,t,\rho)\in\mathbb{T}^{d+1}_{s_{m}}\times B(r_{m})$, the
proof is the same as that in the proof of Lemma 2.1. This completes the proof
of the Iterative Lemma. ∎

##  4 Proof of the extended Kolmogorov’s Theorem

Let
$\Psi(x,I)=\lim_{m\rightarrow\infty}\Psi_{1}\circ\cdots\circ\Psi_{m}(x,I),(x,I)\in\mathbb{T}^{d+1}_{s_{0}/2}\times
B(r_{0}/2).$ Recall that in the Iterative Lemma, we have proved:

  * (i)

$\mathbb{T}_{s_{0}}^{d+1}\times
B(r_{0})\supset\cdots\supset\mathbb{T}_{s_{j}}^{d+1}\times
B(r_{j})\supset\cdots\supset\mathbb{T}_{s_{0}/2}^{d+1}\times B(r_{0}/2),$

  * (ii)

$\Psi_{j}:\mathbb{T}^{d+1}_{s_{j}}\times
B(r_{j})\rightarrow\mathbb{T}^{d+1}_{s_{j-1}}\times B(r_{j-1})$ is real
analytic.

Then
$\Psi^{(m)}:=\Psi_{1}\circ\cdots\circ\Psi_{m}:\mathbb{T}_{s_{0}/2}^{d+1}\times
B(r_{0}/2)\subset\mathbb{T}_{s_{m}}^{d+1}\times
B(r_{m})\rightarrow\mathbb{T}_{s_{0}}^{d+1}\times B(r_{0})$ is well-defined
and analytic in $\mathbb{T}_{s_{m}}^{d+1}\times B(r_{m}).$ Also by (3.34) and
(3.35) in iterative lemma, we have

$\parallel(\partial\Psi_{j}(x,I)-Id)(z)\parallel_{s_{j},r_{j}}\leq\varepsilon^{-4b^{\prime}}\epsilon_{j-1}\parallel
z\parallel,$ (4.1)

for $\forall z\in T_{(x,I)}(\mathbb{T}_{s_{j}}^{d+1}\times
B(r_{j}))=\mathbb{C}^{d+1}\times\mathbb{C}^{d},$ where $\partial\Psi_{j}$ is
the tangent map of $\Psi_{j},$ $T_{(x,I)}(\mathbb{T}_{s_{j}}^{d+1}\times
B(r_{j}))$ is the tangent space at point
$(x,I)\in\mathbb{T}_{s_{j}}^{d+1}\times B(r_{j}).$ Observe that

$\partial\Psi^{(m)}=(\partial\Psi_{1}\circ\Psi_{2}\circ\cdots\circ\Psi_{m})(\partial\Psi_{2}\circ\Psi_{3}\circ\cdots\circ\Psi_{m})\cdots(\partial\Psi_{m}).$
(4.2)

Then by (4.1),

$\parallel\partial\Psi^{(m)}(x,I)z\parallel_{s_{m},r_{m}}\leq\Pi_{j=1}^{m}(1+\epsilon_{j-1}\varepsilon^{-4b^{\prime}})\parallel
z\parallel\leq\Pi_{j=1}^{m}(1+\frac{1}{2^{j}})\parallel z\parallel\leq
2\parallel z\parallel,\forall z\in\mathbb{C}^{d+1}\times\mathbb{C}^{d}.$ (4.3)

Thus, for $w=(x,I)\in\mathbb{T}^{d+1}_{s_{0}/2}\times B(r_{0}/2),$

$\displaystyle|\Psi^{(m+1)}(w)-\Psi^{(m)}(w)|$ $\displaystyle=$
$\displaystyle|\Psi^{(m)}(\Psi_{m+1}(w))-\Psi^{(m)}(w)|$ $\displaystyle\leq$
$\displaystyle\parallel\partial\Psi^{(m)}(\Psi_{m+1}(w))\parallel_{s_{0}/2,r_{0}/2}\sup_{w\in\mathbb{T}^{d+1}_{s_{0}/2}\times
B(r_{0}/2)}\parallel\Psi_{m+1}(w)-w\parallel$ $\displaystyle\leq$
$\displaystyle 2\varepsilon^{-4b^{\prime}}\epsilon_{m}\leq\epsilon_{m}^{1/2},$

where we have used (3.8), (3.29) and (4.3) in (4). Write
$\Psi^{(m)}=\Psi^{(1)}+\sum_{l=2}^{m}(\Psi^{(l)}-\Psi^{(l-1)}),$ here
$\Psi^{(1)}:=\Psi_{1}.$ Then

$\displaystyle\parallel\Psi(w)-w\parallel_{s_{0}/2,r_{0}/2}$
$\displaystyle\leq$
$\displaystyle\parallel\Psi_{1}(w)-w\parallel_{s_{0}/2,r_{0}/2}+\sum_{l=2}^{\infty}\parallel\Psi^{(l)}(w)-\Psi^{(l-1)}(w)\parallel_{s_{0}/2,r_{0}/2}\leq\sum_{l=0}^{\infty}\epsilon_{l}^{1/2}\leq\varepsilon.$

This completes the proof of theorem.

##  5 Application to coupled Duffing oscillators

Consider the coupled Duffing oscillators

$\ddot{x_{j}}+x_{j}^{2n+1}+\frac{\partial F(x,t)}{\partial
x_{j}}=0,\;\;j=1,2,\cdots,m,x=(x_{1},...,x_{m})$ (5.1)

where $m,n>0$ are fixed integers, and

$F(x,t)=\sum_{0\leq|\alpha|\leq
2n+1}P_{\alpha}(t)x^{\alpha}=\sum_{0\leq\alpha_{1}+\cdots+\alpha_{m}\leq
2n+1}P_{\alpha_{1}\cdots\alpha_{m}}(t)x_{1}^{\alpha_{1}}\cdots
x_{m}^{\alpha_{m}},\alpha_{j}\in\mathbb{Z}_{+}.$ (5.2)

and $P_{\alpha}(t):\mathbb{T}\rightarrow\mathbb{R}$ is analytic (in this
section, $\mathbb{T}:=\mathbb{R}/\mathbb{Z}$).

Let $\widetilde{A}$ be a large constant. Replacing $x_{j}$ by
$\widetilde{A}x_{j}$ in (5.1), we get

$\widetilde{A}\ddot{x}_{j}+\widetilde{A}^{2n+1}x_{j}^{2n+1}+\widetilde{A}^{-1}\frac{\partial
F(\widetilde{A}x,t)}{\partial x_{j}}=0,$

that is,

$\ddot{x_{j}}+\widetilde{A}^{2n}x_{j}^{2n+1}+\widetilde{A}^{-2}\frac{\partial
F(\widetilde{A}x,t)}{\partial x_{j}}=0,j=1,2,\cdots,m.$ (5.3)

Let

$y_{j}=\widetilde{A}^{-n}\dot{x}_{j},j=1,2,\cdots,m.$

Then

$\dot{y_{j}}=\widetilde{A}^{-n}\ddot{x_{j}}=-\widetilde{A}^{n}x^{2n+1}_{j}-\widetilde{A}^{-(n+2)}\frac{\partial
F(\widetilde{A}x,t)}{\partial x_{j}}.$

Thus, (5.1) can be written as a Hamiltonian system

$\dot{x_{j}}=\frac{\partial H}{\partial y_{j}},\dot{y_{j}}=-\frac{\partial
H}{\partial x_{j}},j=1,2,\cdots,m,$ (5.4)

where

$H=\widetilde{A}^{n}\sum_{j=1}^{m}\left(\frac{1}{2}y_{j}^{2}+\frac{1}{2(n+1)}x_{j}^{2n+2}\right)+\widetilde{A}^{-(n+2)}F(\widetilde{A}x,t).$
(5.5)

Consider an auxiliary Hamiltonian system which is autonomous:

$\dot{x}=\frac{\partial H_{0}}{\partial y},\dot{y}=-\frac{\partial
H_{0}}{\partial
x},H_{0}=\frac{1}{2}y^{2}+\frac{1}{2(n+1)}x^{2n+2},(x,y)\in\mathbb{R}^{2}.$
(5.6)

Let $(x,y)=(u_{0}(t),v_{0}(t))$ be the solution to (5.6) with initial values
$(u_{0}(0),v_{0}(0))=(1,0).$ Then this solution $(u_{0},v_{0})$ is clearly
periodic. Let $T_{0}$ be its minimal period. By energy conservation, we have

$(n+1)v^{2}_{0}(t)+u_{0}^{2n+2}(t)=1,\forall t\in\mathbb{R}.$ (5.7)

Let

$\left\\{\begin{array}[]{ll}x=c^{\alpha}I^{\alpha}u_{0}(\theta T_{0}),\\\
y=c^{\beta}I^{\beta}v_{0}(\theta
T_{0}),\end{array}\right.(\theta,I)\in\mathbb{T}\times\mathbb{R}_{+},$

where $\alpha=\frac{1}{n+2},\beta=1-\alpha,c=\frac{1}{\alpha T_{0}}.$ It is
easy to check $\det\frac{\partial(x,y)}{\partial(\theta,I)}=1.$ That is,
$dx\wedge dy=d\theta\wedge dI.$ Thus

$\Phi:\left\\{\begin{array}[]{ll}x_{j}=c^{\alpha}I^{\alpha}_{j}u_{0}(\theta_{j}T_{0}),\\\
y_{j}=c^{\beta}I^{\beta}_{j}v_{0}(\theta_{j}T_{0}),\end{array}\right.j=1,2,\cdots,m.$
(5.8)

is well-defined and symplectic transformation
$\Phi:(\theta,I)\in\mathbb{T}^{m}\times\mathbb{R}_{+}^{m}\rightarrow\mathbb{R}^{m}\times\mathbb{R}^{m}.$
Moreover,

$H^{(1)}:=H\circ\Phi=\widetilde{H}_{0}(I)+\widetilde{R}(\theta,t,I),$

where

$\widetilde{H}_{0}(I)=\frac{c^{2\beta}}{2(n+1)}\widetilde{A}^{n}\sum_{j=1}^{m}I_{j}^{\frac{2(n+1)}{n+2}},$
$\widetilde{R}(\theta,t,I)=\widetilde{A}^{-(n+2)}F(\widetilde{A}c^{\alpha}I_{1}^{\alpha}u_{0}(\theta_{1}T_{0}),\cdots,\widetilde{A}c^{\alpha}I_{m}^{\alpha}u_{0}(\theta_{m}T_{0}),t).$

Note that $F$ is a polynomial in $x\in\mathbb{R}^{m}$ with degree $2n+1.$ So
$\widetilde{R}(\theta,t,I)$ is a polynomial of $\widetilde{A}$ of degree $n-1$
with its coefficients analytically depending on
$(\theta,t,I)\in\mathbb{T}^{m+1}\times[1,2].$ Let
$\varepsilon=\widetilde{A}^{-1},a=n,b=n-1,$ and

$H_{0}(I)=\frac{c^{2\beta}}{2(n+1)}\sum_{j=1}^{m}I_{j}^{\frac{2(n+1)}{n+2}},R(\theta,t,I)=\varepsilon^{b}\widetilde{R}(\theta,t,I).$

Then

$H^{(1)}=\varepsilon^{-a}H_{0}(I)+\varepsilon^{-b}R(\theta,t,I).$

where $R(\theta,t,I)$ depends on $\varepsilon$. Note that it is harmless not
to write explicitly the dependence of $R$ on $\varepsilon$. Applying Theorem
1.1 and Remark 3 to $H^{(1)}$, we have that there exists a subset
$\tilde{J}_{0}(\widetilde{A})\subset[1,2]^{m}$ with $\text{Leb
}\,\tilde{J}_{0}(\widetilde{A})\geq 1-(\log\widetilde{A})^{-c_{0}}$ (some
$c_{0}>0$) such that for any solution to (5.4) starting from
$I(0)=(I_{1}(0),\cdots,I_{m}(0))\in\tilde{J}_{0}(\widetilde{A})$ is quasi-
periodic with rotational frequency $(\omega(I(0)),1)$, where
$\omega=(\omega_{j}:\,j=1,...,m)$ and
$\omega_{j}=\varepsilon^{-a}\frac{\partial H_{0}(I(0))}{\partial I_{j}}.$
Returning to the coordinates $(x_{i},\dot{x}_{i}:i=1,...,m)$, we have that for
any large $\widetilde{A}$ there exist sets $\Theta_{\widetilde{A}}$ and
$\tilde{\Theta}_{\widetilde{A}}$ with
$\Theta_{\widetilde{A}}\subset\tilde{\Theta}_{\widetilde{A}}\subset[C_{1}\widetilde{A},C_{2}\widetilde{A}]^{2m}$
and
$\text{Leb}\,\tilde{\Theta}_{\widetilde{A}}=1,\;\text{Leb}\,\Theta_{\widetilde{A}}\geq
1-\left(\log\widetilde{A}\right)^{-C_{0}}$ such that for any solution to (5.1)
starting from $(x_{j}(0),\dot{x}_{j}(0):j=1,...,m)\in\Theta_{\widetilde{A}}$
is quasi-periodic with rotational frequency $(\omega(I(0)),1)$. Note that any
quasi-periodic solution is bounded. Also note that

$\lim_{\widetilde{A}\to\infty}\frac{\tilde{\Theta}_{\widetilde{A}}}{\Theta_{\widetilde{A}}}=1$

Thus, we see that (5.1) is almost Lagrangian stability. This completes the
proof of Theorem 1.2.

## Acknowledgements

The work was supported in part by National Natural Science Foundation of China
(Nos. 11771093, 12071254).

## References

  * (1) G. Duffing, Erzwungene schwingungen bei veranderlicher eigenfrequenz, F. Vieweg und Sohn, Braunschweig, 1918.
  * (2) F. C. Moon, P. J. Holmes, A magnetoelastic strange attractor. J. Sound Vib., 65 (2), 285-296, 1979.
  * (3) F. C. Moon, P. J. Holmes, Addendum: a magnetoelastic strange attractor. J. Sound Vib., 69 (2), 339, 1980.
  * (4) G. Morris, A case of boundedness of Littlewood s problem on oscillatory differential equations, Bull. Austral. Math. Soc., 14, 71-93, 1976.
  * (5) R. Dieckerhoff, E. Zehnder, Boundedness of solutions via twist theorem, Ann. Scula. Norm. Sup. Pisa, 14, 79-95, 1987.
  * (6) J. Moser, On invariant curves of aera-preserving mapping of annulus, Nachr. Akad. Wiss. Gottingen Math.Phys., 2, 1-20, 1962.
  * (7) M. Levi, Quasi-periodic motions in superquadratic time-periodic potentials, Comm. Math. Phys., 143, 43-83, 1991.
  * (8) B. Liu, Boundedness for solutions of nonlinear Hill’s equations with periodic forcing terms via Moser’s twist theorem, J. Diff. Eqs, 79, 304-315, 1989.
  * (9) B. Liu, Boundedness for solutions of nonlinear periodic differential equations via Moser’s twist theorem, Acta Math Sinica (N.S.), 8, 91-98, 1992.
  * (10) X. Yuan, Invariant Tori of Duffing-type Equations, Advances in Math (China), 24, 375-376, 1995.
  * (11) X. Yuan, Invariant Tori of Duffing-type Equations, J. Diff. Equations, 142(2), 231-262, 1998.
  * (12) H. Rüssmann, On the existence of invariant curves of twist mappings of an annulus, Geometric Dynamics. Lecture Notes in Mathematics, 677-718, 1983.
  * (13) R. Jothimurugan, K. Thamilmaran, S. Rajasekar, M. A. F. Sanjuán, Multiple resonance and anti-resonance in coupled Duffing oscillators, Nonlinear Dyn., 83, 1803-1814, 2016.
  * (14) István Z. Kiss, Resonance clustering in globally coupled electrochemical oscillators with external forcing, Phy. Rev. E, 77, 046204, 2008.
  * (15) Agnessa Kovaleva, Capture into resonance of coupled Duffing oscillators, Phy. Rev. E, 92, 022909, 2015.
  * (16) M. G. Clerc, S. Coulibaly, M. A. Ferré, R. G. Rojas, Chimera states in a Duffing oscillators chain coupled to nearest neighbors, Chaos, 28, 083126, 2018.
  * (17) Prasun Sarkar, Deb Shankar Ray, Vibrational antiresonance in nonlinear coupled systems, Phy. Rev. E, 99, 052221, 2019.
  * (18) Amit Kumar Chatterjee, Anupam Kundu, Manas Kulkarni, Spatiotemporal spread of perturbations in a driven dissipative Duffing chain: An out-of-time-ordered correlator approach, Phy. Rev. E, 102, 052103, 2020.
  * (19) J. Shena, N. Lazarides, J. Hizanidis, Multi-branched resonances, chaos through quasiperiodicity, and asymmetric states in a superconducting dimer, Chaos 30, 123127, 2020.
  * (20) J. P. Deka, A. K. Sarma, A. Govindarajan, M. Kulkarni, Multifaceted nonlinear dynamics in PT-symmetric coupled Liénard oscillators, Nonlinear Dyn., 100, 1629-1640, 2020.