\documentstyle[12pt,pstricks,amsfonts,epsfig,graphicx,subfigure]{article}

%\usepackage[textures]{epsfig}

\setlength{\hoffset}{-1.cm}

\setlength{\voffset}{-0.5cm}

\setlength{\textwidth}{17cm}

\setlength{\textheight}{20cm}

\setlength{\parindent}{8mm}

\setlength{\footskip}{2.truecm}
\frenchspacing

\newcommand\beq[1]{ \begin{equation}\label{#1} }

\newcommand{\eeq}{ \end{equation} }

\newcommand\equ[1]{{\rm (\ref{#1})}}

\def\norm#1{\left \vert#1\right \vert}
\def\Norm#1{\Vert #1\Vert }
\def \R  {\bf R}
\def \T  {\bf T}
\def \Z  {\bf Z}
\def \N  {\bf N}
\def \C  {\bf C}




\begin{document}





\title{\bf From regular to chaotic motions through the work of Kolmogorov}





\author{Alessandra Celletti$^{(1)}$, Claude Froeschl\'e$^{(2)}$ and Elena Lega$^{(3)}$}



\maketitle



{$^{(1)}$ Dipartimento di Matematica,
Universit\`{a} di Roma ``Tor Vergata'',
Via della Ricerca Scientifica, 00133 Roma, Italy,
e--mail: celletti@mat.uniroma2.it}


{$^{(2)}$ Observatoire de Nice, 

BP 229, F--06304 Nice Cedex 4 (France), 

e--mail: claude@obs-nice.fr}



{$^{(3)}$ Observatoire de Nice, 

BP 229, F--06304 Nice Cedex 4 (France), 

e--mail: elena@obs-nice.fr}



\section{Introduction}
We know that the solar system survived catastrophic events since 5 billion
years,
which corresponds to the estimated age of our star. Although the Sun will
become a red
giant in about 5 billion years and it will destroy the regularity of the
planetary motions,
the question of the stability of the solar system within this time interval
remains open.

This problem has been faced since antiquity and it motivated the
development of many studies in mathematical physics. In this context, the
so--called
\sl perturbation theory \rm was developed by Laplace, Leverrier, Delaunay
and Lagrange,
with the aim to investigate the dynamics of the planetary motion. To be
more precise,
let us start with the simplest model of planetary evolution, i.e. the
two--body
problem, consisting in the study of the motion of a planet under the
gravitational
influence of the Sun. As is well known, the solution is provided by
Kepler's laws,
which guarantee that (in this approximation) the planet orbits the Sun on
an
elliptic orbit, the Sun being located at one of the foci. The two--body
problem
is the typical example of an \sl integrable \rm system, since its solution
can be
explicitly provided. In order to get a more realistic model, we should
include
the gravitational attraction of the other bodies of the solar system.
Therefore,
let us start by considering the \sl three--body problem, \rm where the
effect of
a third body is included. Despite many efforts performed by a lot of
scientists in
the XIX century, a mathematical solution of this problem cannot be found:
indeed,
H. Poincar\'e proved the non--integrability of the three-body
problem (\cite{Po}). However, if we consider that the action of the third
body
(whatever it is: planet, satellite or comet) is small compared to the
force exerted by
the Sun, we can include the three--body problem in the class of the \sl
nearly--integrable \rm
systems. In our example, the integrable approximation is provided by the
two--body problem and
the strength of the nonlinearity is relatively small (being proportional
to the mass ratio of
the primaries, i.e. the third body and the Sun). In this framework,
perturbation theory provides
useful tools to give (in particular cases) an approximate solution of the
equations of motion.
However, such theory does not allow to draw definite conclusions about the
stability of the dynamics,
since the mathematical series involved in the construction of the
approximate solution are generally
divergent.

It was only in the '50s that A.N. Kolmogorov made a major breakthrough in
the study of nearly--integrable systems.
He proved (\cite{K}) that for slightly perturbed systems, under quite
general assumptions,
some regions of the phase space remain regular. We stress that
Kolmogorov's theorem
can be used in a constructive way to prove the stability of
motions. A few years later, the pioneering
work of Kolmogorov was extended by V.I. Arnold (\cite{Ar1}) and J. Moser
(\cite{M}),
giving rise to the so-called KAM theory. Though the original estimates of
the
KAM theorem do not provide rigorous results in agreement with the
physical predictions, recent implementations of computer-assisted
proofs show that KAM theory can be
efficiently used in model problems, yielding results in agreement with the
experimental
observations.  KAM theory was the starting point for a broad investigation
of the
stability of nearly-integrable systems, which is briefly reviewed in \S2.
Moreover, it opened a new field of research, related to the development of
numerical experiments.

In this framework, the simulations performed by M. H\'enon and C. Heiles
(\cite{HenHeil64})
filled the gap between theory and experiments, opening a bridge toward the
understanding of the global study of the dynamical properties of
the whole phase space.
In  particular, the concept of chaos makes its appearance in a wide range
of
physical problems. The extent of chaotic motions is provided by the
computation of Lyapunov's exponents, which allow to measure the
divergence of nearby trajectories. This concept was recently refined to
investigate the behaviour of weakly chaotic motions, through the
implementation of frequency analysis and Fast Lyapunov Indicators
(hereafter, FLI).
In \S3 we review some applications of such methods to investigate the
topology
of the phase space. Moreover, the computation of the FLIs allowed to
study the transition from Chirikov's to Nekhoroshev's regime and therefore
it provides results about diffusion in Hamiltonian systems.



\section{Stable motions}
\subsection{Integrable and non--integrable systems}

When dealing with real physical problems, it is rather exceptional to find a 
mathematical solution which describes the dynamics of the system. If the solution 
cannot be found, the typical approach is to diminish the complexity of the problem 
by introducing a reasonable set of assumptions. If a solution of the most simplified model 
can be found, then one can proceed to reintroduce the ingredients which make the 
model closer to reality. However, such procedure can fail even at the basic step, when the 
simplest model is considered. 

A familiar example of this approach is provided by the investigation 
of the stability of the solar system. The reality consists of 
9 planets orbiting around the Sun, a lot of satellites moving about the planets, thousands of small
objects like comets and asteroids. Moreover, one should take into account the oblateness of the 
solar system bodies, the existence of rings around the major planets, tides, solar wind effects
and so on. Without proceeding in the discouraging list of events which contribute to the
description of the real world, we start our investigation by reducing the problem to the study of
the simplest model arising in Celestial Mechanics: the two--body problem. Consider,
for example, the motion of an asteroid about the Sun, neglecting all additional effects like 
the gravitational influence of the other planets. The solution of this simplified problem dates
back to J. Kepler, according to which the asteroid revolves around the Sun on an elliptic orbit  
and the motion can be described by elementary formulae. The two--body problem is the archetype
of an \sl integrable \rm system, for which an explicit solution can be found. 

According to our strategy we increase the difficulty of the problem, just adding the effect of a
third body, for example Jupiter, the largest planet of the solar system. The three--body problem
has been investigated since the XVIII century by several astronomers and mathematicians. However,
an explicit solution cannot be found. Indeed, H. Poincar\'e (\cite{Po}) proved that the three--body problem
is
\sl non--integrable, \rm since there are not sufficiently many independent integrals (namely, quantities
which are constant during the evolution of the dynamics), allowing to derive mathematical
expressions which describe the motion of an asteroid under the Newtonian attraction of the Sun and
Jupiter. Since the mass of Jupiter is much smaller than that of the Sun (the mass ratio
amounts to about $10^{-3}$), the asteroid--Sun--Jupiter problem is close to the asteroid--Sun integrable
system whose solution is provided by Kepler's laws. Therefore the effect of Jupiter can be
considered as a small perturbation and the three--body problem enters the class of the so--called 
\sl nearly--integrable \rm systems. In this case, even if the complete answer cannot be given, one
can proceed trying to find an \sl approximate \rm solution of the equations of motion, which
describe the dynamics with \sl sufficient \rm accuracy. The development of this kind of approach is known as
\sl perturbation theory, \rm which we shall describe in the next section. 

\vskip.1in 

In order to make these concepts quantitative, let us consider a mechanical system 
with $n$ degrees of freedom, described by the Hamiltonian function 
$\tilde H({\underline p},{\underline q})$, where ${\underline p}\in{\bf R}^n$, ${\underline
q}\in{\bf R}^n$. If the system is integrable, there exists a canonical change of variables 
${\cal C}:({\underline I},{\underline \varphi})\in{\bf R}^n\times{\bf T}^n \rightarrow ({\underline
p},{\underline q})\in{\bf R}^{2n}$ (with ${\bf T}\equiv {\bf R}/2\pi{\bf Z}$), such that  the transformed
Hamiltonian becomes 
$$
H\circ{\cal C}({\underline I},{\underline \varphi})=h({\underline I})\ , 
$$
having the property to depend only on the variables ${\underline I}$. The coordinates $({\underline
I},{\underline
\varphi})$ are known in the literature as \sl action--angle \rm variables (\cite{Arbook}). The
integrability of the system can be immediately recognized by writing down Hamilton's equations. In fact,
denoting by 
$$
{\underline\omega}={\underline\omega}({\underline I})=
{{\partial h({\underline I})}\over {\partial{\underline I}}}
$$ 
the \sl frequency \rm or \sl rotation
number \rm of the system, one has  
\begin{eqnarray*}
\dot {\underline I}&=&{{\partial h({\underline I})}\over {\partial{\underline
\varphi}}}={\underline 0}\nonumber\\
\dot {\underline \varphi}&=&{{\partial h({\underline I})}\over {\partial{\underline I}}}={\underline\omega}({\underline I})\ . \nonumber\\
\end{eqnarray*}
Therefore, the vector ${\underline I}$ is constant during the motion, say 
${\underline I}={\underline I}_0$, while from the second equation we have ${\underline
\varphi}={\underline\omega} ({\underline I}_0)t+{\underline \varphi}_0$, where $({\underline I}_0,{\underline \varphi}_0)$ denote the initial conditions. 

A nearly--integrable system can be conveniently described in terms of a Hamiltonian of the form 
\begin{equation}
H({\underline I},{\underline \varphi})\ =\ h({\underline I})+\varepsilon f({\underline I},{\underline
\varphi})\ , 
\label{H}
\end{equation}
where $h({\underline I})$ represents the integrable part, while $\varepsilon f({\underline I},{\underline\varphi})$ is the  perturbation whose size, measured by the parameter 
$\varepsilon$, is supposed to be small. 

The dynamics associated to the three--body problem can be described in terms of a nearly--integrable
Hamiltonian system of the previous form. Indeed, referring to the example of the motion of the asteroid under
the attraction of the Sun and Jupiter, the integrable part represents the two--body asteroid--Sun
interaction  (whose solution is given by Kepler's laws), while the perturbation is due to the gravitational
influence of Jupiter; in particular, the perturbing parameter $\varepsilon$ represents the Jupiter-Sun mass
ratio.  Though an explicit solution of this problem cannot be found, perturbation theory can provide useful
tools to predict the motion within good accuracy. 



\subsection{Perturbation theory} 
Celestial Mechanics stimulated the development of perturbation theory. To give an example, thanks to
unexplained perturbations on the motion of Uranus, the
discovery of the planet Neptune was anticipated theoretically by J. Adams and U. Leverrier; their
mathematical computations, performed using perturbation theory,  predicted the position of Neptune with
astonishing accuracy. 

In order to explore in more detail the fundamentals of perturbation theory (\cite{Arbook}, \cite{Arenc},
\cite{G}, \cite{SM}), let us start by considering a Hamiltonian function of the form \equ{H}. Neglecting the
perturbation, we can solve explicitly the equations corresponding to the Hamiltonian $h({\underline I})$; this provides an $\varepsilon$--approximation of the true solution up to 
a time of the order of $1/\varepsilon$. The goal of
perturbation theory is to introduce a suitable change of coordinates such that we can integrate the
equations of motion for  a longer time. In particular, we describe a simple algorithm which allows to remove
the perturbation to  order $\varepsilon^2$. More precisely, we define a convenient change of variables, say 
${\cal C}:({\underline I'},{\underline \varphi'})\rightarrow ({\underline I},{\underline \varphi})$ with
the property that the transformed Hamiltonian $H'({\underline I'},{\underline \varphi'})$ takes the form 
\begin{equation}
H'({\underline I'},{\underline \varphi'})\ =\ H\circ{\cal C}({\underline I'},{\underline \varphi'})\ \equiv \ 
h'({\underline I'})+\varepsilon^2 f'({\underline I'},{\underline\varphi'})\ , 
\label{H'}
\end{equation}
with suitable functions $h'$ and $f'$. We remark that the perturbing function has now been removed to orders 
$\varepsilon^2$ and therefore Hamilton's equations associated to \equ{H'} can be integrated up to a time
of order $1/\varepsilon^2$. Since $\varepsilon$ is a small quantity ($0<\varepsilon<<1$) we get the solution
of the motion for a longer time. This idea was developed since the XVIII century to compute the dynamics
of the bodies of the solar system with more precision than Kepler's laws. The ingredients needed to implement
perturbation theory are very elementary and are based on the construction of a coordinate's
transformation, a Taylor series development as a function of the perturbing parameter $\varepsilon$ and a
Fourier series expansion to determine explicitly the change of coordinates. 

It is rather instructive to derive a first order perturbation theory. We ask the reader to make a little
effort in trying to  follow the mathematical formulae, so that she/he can control the basics of the
computations, which allow to make refined predictions of the solar system dynamics. Let us start by
introducing the canonical change  of variables $({\underline I},{\underline \varphi})\rightarrow
({\underline I'},{\underline \varphi'})$ by means of the system of equations 
\begin{eqnarray*}
{\underline I}&=&{\underline I}'+\varepsilon 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline \varphi}}}\nonumber\\
{\underline \varphi}'&=&{\underline \varphi}+\varepsilon 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline I}'}}\ , \nonumber\\
\end{eqnarray*}
where the \sl unknown \rm function $\Phi({\underline I}',{\underline \varphi})$ is usually denoted as the \sl
generating function. \rm We remark that $\Phi$ depends on the old angle variables ${\underline \varphi}$ and
on the new action variables ${\underline I}'$, as usually defined in textbooks of classical mechanics (see,
e.g., \cite{Arbook}). In order to give explicit formulae for the transformed Hamiltonian, let us start by 
splitting the perturbing function $f({\underline I},{\underline\varphi})$ into a part 
$f_0({\underline I})={1\over {(2\pi)^n}}\int_{{\bf T}^n}f({\underline I},{\underline \varphi})d{\underline
\varphi}$, depending only on the variables ${\underline I}$ and a remainder function $\tilde f({\underline I},  {\underline \varphi})$ defined by the expression $\tilde f({\underline I},{\underline \varphi})\equiv 
f({\underline I},{\underline\varphi})-f_0({\underline I})$.  
Next, we insert the change of variables in the Hamiltonian \equ{H} and we expand in Taylor series 
up to the second order around $\varepsilon=0$; more precisely, we obtain 
\begin{eqnarray*}
&&h({\underline I}'+\varepsilon 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline \varphi}}})+
\varepsilon f({\underline I}'+\varepsilon 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline \varphi}}},\varphi)
\nonumber\\
&=&h({\underline I}')+{\underline\omega}({\underline I}')\cdot \varepsilon 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline \varphi}}}+
\varepsilon f_0({\underline I}')+\varepsilon \tilde f({\underline I}',{\underline
\varphi})+O(\varepsilon^2)\ , \nonumber\\
\end{eqnarray*}
where we recall that ${\underline\omega}\equiv {{\partial h}\over {\partial {\underline I}}}$. Since we want
that the new Hamiltonian is integrable up the second power of $\varepsilon$, we need to kill the first order
term in $\varepsilon$ by requiring that the function $\Phi$ satisfies the equation 
\begin{equation}
{\underline\omega}({\underline I}')\cdot 
{{\partial \Phi({\underline I}',{\underline \varphi})}\over {\partial{\underline \varphi}}}+
\tilde f({\underline I}',{\underline \varphi})={\underline 0}\ . 
\label{phi}
\end{equation}
The above equation is well--posed, since its average over the angle variables is zero. 
If we succeed to find a solution of \equ{phi}, then we immediately recognize that the new integrable
Hamiltonian can be written as 
$$
h'({\underline I}')\ =\ h({\underline I}')+\varepsilon f_0({\underline I}')\ , 
$$
whose associated Hamilton's equations provide the solution of the motion up to $O(\varepsilon^2)$. 
The equation \equ{phi} can be used to derive an explicit expression for the generating function. More
precisely, let us expand the functions $\Phi$ and $\tilde f$ in Fourier series as 
\begin{eqnarray*} 
\Phi({\underline I}',{\underline \varphi})&=&\sum_{{\underline m}\in{\bf Z}^n\backslash\{{\underline 0}\}} 
\hat\Phi_{{\underline m}}({\underline I}')\ e^{i{\underline m}\cdot {\underline \varphi}}\ ,\nonumber\\
\tilde f({\underline I}',{\underline \varphi})&=&\sum_{{\underline m}\in{\bf Z}^n\backslash\{{\underline
0}\}} \hat f_{{\underline m}}({\underline I}')\ e^{i{\underline m}\cdot {\underline \varphi}}\ .\nonumber\\
\end{eqnarray*}
Inserting the Fourier expansion inside \equ{phi} we obtain 
$$
i\ \sum_{{\underline m}\in{\bf Z}^n\backslash\{{\underline 0}\}}  {\underline\omega}({\underline I}')\cdot 
{\underline m}\ \hat\Phi_{{\underline m}}({\underline I}')\ e^{i{\underline m}\cdot {\underline \varphi}}
=-\sum_{{\underline m}\in{\bf Z}^n\backslash\{{\underline 0}\}} \hat f_{{\underline m}}({\underline I}')\
e^{i{\underline m}\cdot {\underline \varphi}}\ , 
$$
from which we derive 
$$
\hat\Phi_{{\underline m}}({\underline I}')=-{{\hat f_{{\underline m}}({\underline I}')}\over 
{i\ {\underline\omega}({\underline I}')\cdot {\underline m}}}\ . 
$$
Casting together the above formulae, we obtain that the generating function is given by the expression 
$$
\Phi({\underline I}',{\underline \varphi})= i\ \sum_{{\underline m}\in{\bf Z}^n\backslash\{{\underline 0}\}} 
{{\hat f_{{\underline m}}({\underline I}')}\over {\ {\underline\omega}({\underline I}')\cdot 
{\underline m}}}\ e^{i{\underline m}\cdot {\underline \varphi}}\ . 
$$
We strongly remark that the algorithm presented above is constructive in the sense that it allows to give
explicit formulae for the computation of the transformed Hamiltonian as well as of the canonical change of
variables. However, there could be a problem in the above implementation whenever the terms
${\underline\omega} ({\underline I}')\cdot {\underline m}$ are zero, meaning that the frequency vector is
rationally dependent. In such case, the generating function is not defined and the algorithm fails. Indeed, 
if the vector ${\underline\omega}$ is rationally independent, the divisors are not identically zero, but 
they can become arbitrarily small and they can lead to the
divergence of the series defining the generating function. The \sl small divisor \rm problem is the obstacle
which prevents the iteration of the  algorithm described above to higher orders in $\varepsilon$. 
If the divisors are very small, the Fourier coefficients defining the generating function become 
very big and it may happen that the series does not converge. However, the theory founded by A.N. 
Kolmogorov in 1954 allows to iterate infinitely many times the above procedure and, 
under suitable assumptions, to prove the existence of \sl stable \rm motions.  


\subsection{The Kolmogorov--Arnold--Moser theorem}
The success of perturbation theory crucially depends on the character of
the  frequency vector ${\underline\omega}$. Let us focus on the case of
${\underline\omega}=(\omega_1,...,\omega_n)$ being rationally independent, namely there do not exist a
vanishing linear combination with rational coefficients  (unless all the coefficients are zero). A \sl
quasi--periodic \rm motion is a solution that can be expressed (as the time
$t$ varies) as a function of the form $F=F(\omega_1 t,...,\omega_n t)$,  where $F(\varphi_1,...,\varphi_n)$
is a multiperiodic function, $2\pi$--periodic in the angles $\varphi_i$ and 
${\underline\omega}$ is rationally independent. 
In order to illustrate the difference between periodic and quasi--periodic motions, let us consider the
dynamics on a two--dimensional torus (i.e., fix $n=2$), which is graphically represented by a life--belt. A
point on the torus is located introducing two angles $\varphi_1$, $\varphi_2$ (with $0\leq \varphi_1\leq
2\pi$, 
$0\leq \varphi_2\leq 2\pi$) measured on the parallel and on the meridian of the torus. If we assume that the
evolution is given by $\varphi_1(t)=\omega_1t+\varphi_1(0)$, $\varphi_2(t)=\omega_2t+\varphi_2(0)$, the
difference between periodic and quasi--periodic motions is related to the rational or irrational
character of the quantity ${\omega_1\over \omega_2}$. In particular, if ${\omega_1\over \omega_2}$ is a
rational number, say ${\omega_1\over \omega_2}={p\over q}$ with $p,q\in{\bf Z}_+$, then the motion is
periodic and it retraces the same steps. On the contrary, if ${\omega_1\over \omega_2}$ is irrational, the
evolution never returns to the initial position and it can be shown that the trajectory is everywhere dense
on the torus, defining a quasi--periodic motion (\cite{Ar2}). 

As described in section 2.2, perturbation theory concerns the fate of the solutions of an integrable system
under a small perturbation which preserves the Hamiltonian structure. Standard theories approached this
problem by fixing the initial conditions and by investigating the consequent evolution. 
Instead of looking for the stability of the solution with
preassigned initial conditions, Kolmogorov changed point of view by investigating the stability of the dynamics
with a fixed frequency ${\underline\omega}$. At the International Congress of Mathematicians held in Amsterdam in
1954, Kolmogorov announced a result which marked a milestone in the development of the stability theory of
nearly--integrable Hamiltonian systems (\cite{K}). 

Kolmogorov's theorem can be stated as follows: 
consider a real--analytic Hamiltonian function of the form \equ{H} and fix a rationally independent frequency
vector
${\underline\omega}$; under suitable assumptions on the unperturbed Hamiltonian $h$ and on the frequency  
${\underline\omega}$, if the strength of the perturbation (namely the parameter $\varepsilon$) is
sufficiently small, there exists an invariant torus on which a quasi--periodic motion with frequency
${\underline\omega}$  takes place. Moreover, Kolmogorov's theorem states that the collection of such
invariant tori form a set of positive measure in the phase space. 
In \cite{K} Kolmogorov gave also an outline of the proof: such scheme turned out to be particularly rich 
and useful. 
In 1963 V.I. Arnold (\cite{Ar1},
\cite{Ar2}) published a detailed alternative proof of Kolmogorov's theorem, while in 1962  J. Moser
(\cite{M}) developed a proof about the existence of invariant curves for smooth enough area--preserving
mappings in the plane. Since then, the overall theory is known with the acronym of KAM theory. 

\vskip.1in 

In order to illustrate KAM theory in more detail, let us consider a Hamiltonian function of the form
\equ{H}.  As far as the integrable system is concerned (namely setting $\varepsilon=0$), we have seen that
Hamilton's equations lead to quasi--periodic motions where the actions are constant, say ${\underline I}={\underline I}_0$  with ${\underline\omega}({\underline I}_0) \equiv {\underline\omega}_0$, 
while the angles vary linearly with the time (${\underline \varphi}(t)={\underline\omega}_0t+{\underline
\varphi}_0$). The assumptions under which KAM theory can be applied are very general. The first hypothesis
concerns  the unperturbed Hamiltonian $h({\underline I})$ which we assume to be \sl non--degenerate, \rm in
the sense that the  determinant of the Hessian matrix is different from zero 
\footnote{Notice that the non--degeneracy condition implies that the
 frequencies satisfy the relation:
$\det {{\partial {\underline\omega} ({\underline I})}\over {\partial {\underline I}}}
\not=0$.}:
\begin{equation}
\det {{\partial^2 h({\underline I})}\over {\partial {\underline I}^2}}\not=0\ ,\qquad \forall\ 
{\underline I}\in {\bf R}^n\ . 
\label{c1}
\end{equation}
The second assumption concerns the frequency vector: in order to avoid the small--divisor problem, one 
requires that ${\underline\omega}_0$ satisfies a strong irrationality condition, namely that 
${\underline\omega}_0$ fulfills the \sl diophantine \rm inequalities: 
\begin{equation}
|{\underline\omega}_0\cdot {\underline m}|^{-1}\ \leq\ C|{\underline m}|^\tau\ , \qquad \forall\ {\underline
m}\in {\bf Z}^n\backslash\{{\underline 0}\}\ , 
\label{c2}
\end{equation}
where $C$, $\tau$ are positive constants. 
We remark that the above assumption is more stringent than requiring a \sl non--resonance \rm condition 
of the form $|{\underline\omega}_0\cdot {\underline m}|\not=0$ for any ${\underline m}\in
{\bf Z}^n\backslash\{{\underline 0}\}$. 

Since the diophantine condition \equ{c2} plays a key role in the statement of the KAM theorem, we 
want to make it clearer through a concrete example. Take $n=2$ and assume that the frequency 
vector has the form ${\underline\omega}_0\equiv (\gamma,1)$ with $\gamma$ being a real number. Let us
introduce the
\sl continued  fraction \rm representation of the number $\gamma$ as the sequence of positive integers
$\{a_0,a_1,a_2,a_3...\}$,  such that 
$$
\gamma=a_0+{1\over {a_1+{1\over {a_2+{1\over {a_3+...}}}}}}\ .
$$
Using a standard notation, we write $\gamma=[a_0;a_1,a_2,...]$. A simple computation shows that if 
$\gamma$ is a rational number, then its continued fraction is composed by a finite number of terms: 
there exists a positive integer $N$ such that $\gamma=[a_0;a_1,a_2,...,a_N]$. 
On the contrary, if $\gamma$ is an irrational number, then the continued fraction representation is infinite. 
Number theory guarantees that the diophantine condition \equ{c2} with $\tau=1$ is satisfied, for example, by
the  so--called \sl noble \rm numbers, i.e. irrational numbers whose continued fraction is definitely
1: there exists a positive integer $M$ such that $\gamma=[a_0;a_1,a_2,...,a_M,1,1,1,...]$. Obviously, the
most "noble" irrational number is represented by a continued fraction composed by all 1's;
this number, known since antiquity, is the \sl golden ratio \rm ${{\sqrt{5}+1}\over 2}=[1;1,1,1,...]$. The
golden ratio satisfies  condition \equ{c2} with the smallest diophantine constant $C$, being
$C={{3+\sqrt{5}}\over 2}$. Finally, let us mention that each irrational number can be approximated by a
sequence of rational numbers $\{{p_k\over q_k}\}_{k\in{\bf Z}}$ provided by the successive truncations of
the continued fraction, i.e. 
$$
{p_1\over q_1}\ =\ a_0+{1\over a_1}\ ,\qquad {p_2\over q_2}\ =\ a_0+{1\over {a_1+{1\over {a_2}}}}\ ,\qquad 
{p_3\over q_3}\ =\ a_0+{1\over {a_1+{1\over {a_2+{1\over {a_3}}}}}}\ ... 
$$
For the golden ratio, the rational approximants are given by the ratio of the so-called Fibonacci's
numbers and the sequence $\{1,{2\over 1},{3\over 2},{5\over 3},{8\over 5},{{13}\over 8},{{21}\over
{13}},...\}$ converges to 
${{\sqrt{5}+1}\over 2}$. 

\vskip.1in 

Coming back to the description of the content of KAM theory, let us consider a Hamiltonian system described
by
\equ{H}, assuming that the conditions \equ{c1} and \equ{c2} are satisfied. We ask ourselves whether there
still   exists for the perturbed problem (namely taking
$\varepsilon\not=0$) an invariant torus on which  a quasi--periodic motion with frequency ${\underline\omega}_0$
takes place. KAM theory yields a positive answer to this question, provided that the size of the perturbing
parameter $\varepsilon$ is sufficiently small,  say $\varepsilon\leq \varepsilon_{KAM}({\underline\omega}_0)$. The
smallness condition is required in   order to guarantee the convergence of the series developments involved in the
proof.  We remark that for low--dimensional iso--energetically non--degenerate 
Hamiltonian systems, the existence of a KAM torus provides a 
strong stability property. More precisely, if $n=2$ the phase space has dimension 4 and the 
constant energy surfaces have dimension 3. Therefore, the 2--dimensional KAM tori foliate the 
constant energy levels providing the stability of the motion, due to the uniqueness of the 
solutions of ordinary differential equations. This property ceases to be valid as $n\geq 3$; 
indeed, already if $n=3$ the dimensions of the phase space and of the constant energy surfaces are,  
respectively, 6 and 5. Thus, the 3--dimensional KAM tori do not provide a separation 
of the phase space into invariant regions as it is the case for the 2--dimensional Hamiltonian 
systems. 

\vskip.1in 

Before having an insight on the proof of the KAM theorem, let us briefly describe the qualitative behaviour 
of the KAM torus as the perturbing parameter is varied. If $\varepsilon$ is very small, the 
invariant torus lives close to the location of the unperturbed torus, but it becomes more and more 
displaced and deformed as the perturbing parameter increases. At a certain value of 
$\varepsilon$, say $\varepsilon=\varepsilon_c({\underline\omega}_0)$, the invariant torus ceases to 
be regular and for $\varepsilon>\varepsilon_c({\underline\omega}_0)$ breaks down. 
In order to estimate the critical
break--down threshold of an invariant torus with frequency ${\underline\omega}_0$, we can use the analytical
approach given by the KAM theorem to get a lower bound $\varepsilon_{KAM}({\underline\omega}_0)$.
Alternatively, we can implement a numerical investigation to determine the break--down value
$\varepsilon_c({\underline\omega}_0)$. Along this direction, the most used numerical method is due to J.
Greene (\cite{Gr}), who provided an efficient algorithm to compute the critical threshold. The original work
(\cite{Gr}) has been developed for a simple discrete system known as the \sl standard mapping, \rm
characterized by a single frequency $\omega_0\in{\bf R}$. Greene's method relies on the conjecture that the
break--down of an invariant surface with frequency
$\omega_0$ is strictly related to the transition from stability to instability of the periodic orbits with
frequency given by the rational approximants $\{{p_k\over q_k}\}_{k\in{\bf Z}}$ converging to $\omega_0$.
The stability character of the periodic orbits can be easily determined by computing the Lyapunov's
exponents; therefore, Greene's method can be easily adapted to numerical computations and it has been
extensively applied to a wide variety of continuous and discrete systems. Partial proofs of Greene's
technique has been presented in \cite{FL}, \cite{MK}. 

\vskip.1in 

Concerning the analytical estimate, 
the idea behind Kolmogorov's theorem is to overcome the small divisor problem by producing a \sl
superconvergent \rm sequence of canonical transformations. We reviewed in section 2.2 the standard perturbation
theory which  allows to transform the initial Hamiltonian, say $H_1=h_1+\varepsilon f_1$, to a Hamiltonian
denoted as $H_2=h_2+\varepsilon^2 f_2$, where the perturbation is, roughly speaking,  of order
$\varepsilon^2$. One can try to iterate  the algorithm to obtain the Hamiltonian $H_3=h_3+\varepsilon^3 f_3$
and, at the
$j$-th step, 
$H_j=h_j+\varepsilon^j f_j$. However, the hindrance due to the convergence of the method cannot be 
in general surmounted. The idea of Kolmogorov consists in the development of a \sl quadratic \rm
method according to which the initial Hamiltonian is transformed to $H_2=h_2+\varepsilon^2 f_2$. 
Iterating one has $H_3=h_3+\varepsilon^4 f_3$ and, at the $j$-th step, $H_j=h_j+\varepsilon^{2^{(j-1)}} f_j$. 
The fact that the perturbation decreases with $\varepsilon$ faster than linearly allows to control 
the small divisors appearing in the sequence of transformations. 

Newton's method for finding the real root of an equation $f(\xi)=0$ offers the opportunity to explain
conveniently the  difference between the linear (standard perturbation theory) and quadratic (KAM
theory) methods (see Figure 1).  Let us start with an initial approximation $\xi_0$ and let
$e_0\equiv|\xi-\xi_0|$ be an estimate  on the initial error. If we proceed to find the solution by a \sl
linear \rm convergence, the  successive approximation $\xi_1$ is determined as the intersection of the
tangent to the curve  at the point $(\xi_0,f(\xi_0))$ with the $\xi$--axis.  Let the derivative of $f$ at
$\xi_0$ be 
$\eta_0\equiv f'(\xi_0)$. Next approximation $\xi_2$ is computed as the abscissa of the line through
$(\xi_1,f(\xi_1))$ with slope $\eta_0$. Iterating this algorithm, the error at step $j$ is
$e_j=|\xi-\xi_j|=O(e_0^{j+1})$, in analogy to standard perturbation theory.  

In the quadratic procedure, the successive approximations are computed as the 
intersection between the $\xi$--axis and the tangent to the function 
computed at the previous step, i.e. 
$$
\xi_{j+1}=\xi_j-{{f(\xi_j)}\over {f'(\xi_j)}}\qquad j=0,1,2...
$$
Let us expand $f(\xi)$ around $\xi_j$; the second order expansion is given by 
$$
0=f(\xi)=f(\xi_j)+f'(\xi_j)(\xi-\xi_j)+{{f''(\xi_j)}\over {2!}}(\xi-\xi_j)^2+O(|\xi-\xi_j|^3)\  , 
$$
which yields 
$$
\xi_{j+1}-\xi={1\over 2!}\ {{f''(\xi_j)}\over {f'(\xi_j)}}(\xi-\xi_j)^2\ . 
$$
Therefore, the error at the step $j$ goes quadratically as  
$$
e_{j+1}=O(e_j^2)=O(e_0^{2^j})\ . 
$$


\begin{figure}[hbt]
  \centering
  \epsfig{file=figure1.eps,width=14cm}
  \caption{Newton's method for finding the root of the equation
$f(\xi)=0$. $a)$ Linear convergence; $b)$ quadratic convergence.}
\end{figure}
The proof of the KAM theorem requires a long set of estimates in order to ensure the convergence of the
method. We refer the reader to \cite{K}, \cite{Ar1}, \cite{M}, \cite{Arenc}, \cite{G}, \cite{SM} for the
detailed proof. 


\subsection{The stability of a model associated to the three--body problem}
The study of the stability of the three--body problem in Celestial Mechanics has a long tradition. A
mathematical approach based on first--order perturbation theory was introduced by Lagrange and Laplace. 
In particular, we focus on the effect of small perturbations on the dynamics of the planets in order to see
whether, after sufficiently long time, collisions or escape to infinity might take place. 
In \cite{Ar2}, the planetary motion in the framework of the three--body and many--body problems was
investigated. Quoting directly V.I. Arnold (\cite{Ar2}, p.125, see also \cite{Ar1}), his result is the
following: \sl "for the majority of the initial conditions under which the instantaneous orbits of the
planets are close to circles lying in a single plane, perturbation of the planets on one another produces,
in the course of an infinite interval of time, little change on these orbits provided the masses of the
planets are sufficiently small". \rm To give concrete estimates, M. H\'enon
(\cite{He})  applied the original version of Arnold's theorem to the three--body problem, allowing to prove 
the existence of invariant tori for values of the perturbing parameter (namely the Jupiter--Sun mass
ratio) less or equal than $10^{-333}$. Such estimate was improved to $10^{-50}$ by applying Moser's theorem.
However, astronomical observations show that the mass of Jupiter compared to that of the Sun is about 
$10^{-3}$. The discrepancy between the KAM estimate and the physical value led to the idea that 
KAM theory is unable to produce effective results. Indeed, at least for some model problems, the refinement
of the KAM estimates  and \sl computer--assisted \rm series expansions provide results which are in
agreement with the  astronomical observations and numerical experiments (\cite{Czamp}, \cite{CCJMP}, \cite{CC0}, 
\cite{CC1}, \cite{CC}). Let
us describe a recent results obtained in
\cite{CC}  concerning the  stability of some asteroids in a 3--body framework (see also \cite{CC1}). We remind that
asteroids are thousands minor bodies of the solar system which populate a belt between the orbits of Mars
and Jupiter. 

Let us introduce a mathematical model, which describes the motion of an asteroid moving under the 
gravitational influence of the Sun and Jupiter. We assume that the three bodies move on the same plane and
that the mass of the asteroid is so small that it does not influence the motion of the primaries. 
Moreover, we assume that Jupiter orbits the Sun on a circular (Keplerian) orbit. This model is known as the
\sl planar, circular, restricted three--body problem. \rm It can be conveniently described in terms of
suitable action--angle coordinates, known as \sl Delaunay's variables \rm (\cite{D}, \cite{Sz}) and the
corresponding Hamiltonian function has two degrees of freedom. The unperturbed frequencies
(namely the frequencies obtained considering the two--body problems provided by the pairs Sun--asteroid or
Sun--Jupiter) are inversely proportional to the periods of revolution about the Sun. We remark that a
resonance condition occurs whenever the ratio of the periods of revolution of the asteroid and of Jupiter
about the Sun is a rational number; there are many examples in the asteroidal belt which meet such
condition, usually denoted as \sl orbit--orbit resonance. \rm 


With reference to section 2.3, since the Hamiltonian system is two--dimensional, we fix a level of energy
(corresponding to the Keplerian motion of the asteroid) in order to have a control on
the dynamics. The stability of the minor body can be obtained by
proving the existence of invariant surfaces which confine the motion of the asteroid on the given energy
level.  Refined KAM estimates for the isoenergetic case have been developed in \cite{CC}.  The proof
involves heavy computations of the series expansion of the solution. To this end, a suitable implementation
on the computer has been performed. We remark that the computer is used so that all numerical errors
introduced by the machine are controlled through the so--called \sl interval arithmetic \rm technique, which
is becoming a widespread tool to prove analytical results with the aid of the computer  (we refer the reader
to \cite{La}, \cite{EW}, \cite{Wi1} and references therein for details about computer--assisted
proofs). 
We want to stress that computer--assisted techniques are a complementary part
in the proof of mathematical results.

The dynamics of the asteroids Iris, Victoria and Renzia have been considered; the application of the 
(isoenergetic, computer--assisted) KAM theory provides their stability for values of the perturbing
parameter less or equal than $10^{-3}$, giving a result in full agreement with the astronomical
prediction. As a conclusion, we stress that this result contributes to show that beside solving the stability
problem from a purely mathematical point of view, KAM theory turns out to be effective also in physical models. 

\section{Unstable motions}
Quasi periodic motions have been investigated in the previous section,
where we showed that the essential features are close to 
the integrable approximation.
Within KAM theory,  nothing is predicted for the initial conditions
which do not satisfy the diophantine inequality \equ{c2}.
In order to have a complete representation of the phase space,
we intend to explore the portion in which a resonance condition of the form  
$\sum_{i=1}^n k_i \omega_i =0$ (with some 
integers $k_1,\ldots ,k_n\in {\bf Z}\backslash \{0\}$) is satisfied in the integrable approximation.
We define the {\sl Arnold's web} as a suitable neighborhood of  resonant
 orbits, whose size depends on the perturbing parameter and on  the order of the resonance, i.e. $\sum_{i=1}^n \vert k_i\vert$.
The trajectories in the Arnold's web can exhibit chaotic features.
 \par
The structure of the Arnold's web was clearly explained already in  
\cite{Arno64}; the numerical investigation of its structure  
has been recently performed on mathematical models 
(\cite{Laskar93}) and on physically interesting systems.  
We remark that the stability analysis in the frame of KAM theory
has touched different fields of physics. To give some examples,
we mention the studies on beam-beam
interactions (\cite{MonthHerrera79}), asteroids diffusion 
(\cite{nesMello97})
 and galactic models (\cite{PapLas98}).
These works have been based on numerical applications using
frequency-map analysis (\cite{Laskar92}).
It is worthwhile to stress that the validity of the
numerical investigations is not only explanatory
or didactic;  indeed, while KAM theorem 
describes the persistence of  regular orbits, a
rigorous proof of the existence of instability or
irregularity in the Arnold web is a delicate problem. \par
In the following sections, we shall  review
a theoretical result due to Nekhoroshev and we will introduce 
 a  numerical tool 
 (\cite{FroLeGo96}),  which allows to  detect and to draw
the Arnold's web in a very detailed way (\cite{FroGuLe00}).
Such analysis allows also to study the 
 transition from a stable to a diffusive regime
(\cite{GuLeFro01}).
 We  conclude by describing a numerical study (\cite{LeGuFro02}) 
of  the so--called Arnold's diffusion (\cite{Arno64}).
\subsection{Nekhoroshev's theorem}
A breakthrough in the solution of the problem of stability of a Hamiltonian
quasi-integrable system is provided by  Nekhoroshev's (\cite{Nek77})
 theorem. Contrary to KAM theorem, Nekhoroshev result extends to an open set of initial conditions; the price to pay is that the 
 stability for infinite times is replaced by the  stability 
over {\sl long} times.
In particular, the  stability time is exponentially long with respect to the
inverse of the perturbing parameter.  Therefore, the concept of long time 
stability is often equivalent in physical problems 
to that of effective stability. To give an example, 
 the infinite time stability of Jupiter has no  physical 
meaning since the solar system will end in some 5 billion years! \par
Nekhoroshev's theorem applies to quasi--integrable Hamiltonians of the
form 
\begin{equation}
H(\underline{I},\underline{ \varphi} )=h(\underline{I})+\varepsilon f(\underline{I},\underline{ \varphi} )\ , \qquad
(\underline{I},\underline{ \varphi}) \in B\times{\bf T}^n \ ,
\label{system}
\end{equation}
where $B$ is an  open subset of ${\R}^n$ and 
${\bf T}^n$ is the $n$--dimensional torus. We assume that $h$ and $f$ 
are analytic functions and that $\varepsilon$ is a small parameter. 
Moreover, the integrable 
approximation $h$ is required to satisfy a suitable geometric condition,
called {\sl steepness} (\cite{Nek77}, see also
\cite{GuMo97},\cite{BeFaGu98}).
In the following, we will assume for simplicity
that $h$ is a  convex function 
(convex functions satisfy the steepness condition).
Under the above hypotheses, the theorem states that, if 
$\varepsilon$ is sufficiently small, any solution 
$(\underline{I}(t),\underline{\varphi} (t))$ of (\ref{system}) satisfies the following exponential 
estimate (\cite{Nek77}, see also
 \cite{BeneGa86}, \cite{Lochak92}, \cite{Poschel93}): 
\begin{equation}
\Norm{\underline{I}(t)- \underline{I}(0)}\leq I _0\, \varepsilon^{a}\qquad  
{\rm for} \quad  |t|\leq t_0\, e^{( \varepsilon_0/\varepsilon )^b}\ ,
\label{Nek}
\end{equation}
where $I_0,t_0,\varepsilon_0,a,b$ are suitable positive constants. In the convex 
case, one can take $a=b=1/(2n)$ (\cite{Lochak92},  \cite{Poschel93}). 
Therefore, for
$\varepsilon < \varepsilon _0$,  the actions  remain close to 
their initial values up to exponentially long times. 
Unfortunately, as often occurs in perturbation theory, 
the purely analytical estimates of the constants 
$I_0,t_0,\varepsilon_0$ are quite unrealistic. This remark motivates
the numerical studies for realistic models of physical interest. \par
The direct numerical check of the instability of the actions
 in general fails, since  actions are
bounded up to an exponentially long time.
 Therefore, alternative  numerical tools have been developed
in the last ten years to investigate the problem in an indirect way
 (\cite{Laskar90}, \cite{Laskar92},  \cite{LeFro96}, \cite{ConVo97}, 
\cite{FroLe97}).
\par
Actually, the hearth of 
 the Nekhoroshev's theorem is that the long time
 stability of the actions occurs within a  specific 
structure of the phase space: the core is made by a large part of KAM tori
while the complementary part is structured by the Arnold's web.
The numerical experiments are based on checking such peculiar structure of
the phase space. 
\subsection{Tools for detecting chaos and order}
In order to explore the phase space, a
careful analysis of a large number of orbits is required.
The classical tool for discriminating between chaotic and ordered
orbits is the computation of the largest Lyapunov exponent, that 
we recall as follows.
Let us consider the differential equations:
\begin{eqnarray}
\left .
\begin{array}{llrr}
 {d \over dt}\underline{ X} \ & = &\underline{ F}(\underline{ X}), & \quad \underline{ X}=(x_1,x_2,...,x_n)
\label{equamotio}
\end{array} \right .
\end{eqnarray}
where $\underline {F}$ are suitable regular functions.
The Lyapunov exponents
are computed by integrating the equations of motion (\ref{equamotio}),
together with  the
corresponding variational
equations:
$$
 {d \underline{ v}\over dt} = ({\partial \underline {F} \over \partial \underline{X}})\, \underline{ v}\ ,
$$
where $\underline{ v}$ is any $n$-dimensional vector.
The largest Lyapunov exponent is defined as  
$$
\gamma = \lim _{t\rightarrow \infty} \ln {||\underline{ v}(t)||\over t} \, ;
$$
let us remark that this limit does not always exist for $v(0)$
belonging to some subspaces of ${\R}^n$ (\cite{Fro84}). 
If (\ref{equamotio}) is of Hamiltonian type, 
 the largest Lyapunov exponent is zero for regular motions and it is 
positive for chaotic orbits;
this property has been largely used to discriminate between chaotic
and ordered motions. However, among regular motions, the Lyapunov 
exponent does not allow to distinguish between KAM tori and resonant 
islands.
When  computing  the Lyapunov exponents, 
 the attention  is focused
on the length of time necessary to get a reliable value
of the limit;
very little importance has been given to the beginning of the integration,
since it was considered a transitory 
regime, depending mainly on the choice
of the initial vector on the  tangent manifold. \par
As  already remarked in  \cite{FroLeGo96}, 
 chaotic, even {\sl weakly} chaotic,
and ordered motion can be discriminated by looking at the value
of $\log ||\underline{v}(T)||$,  where $T$ is relatively small. 
The authors called the above quantity Fast Lyapunov Indicator (hereafter FLI). 
 It turns out  that the FLIs allow also to distinguish 
among regular resonant and non resonant motions 
(see \cite{GuLeFro01} for an analytic explication using a Hamiltonian system).
\par
The precise definition of the FLI can be stated as follows:
  given an initial condition $\underline{ X}(0) \in {\bf R}^{n}$
and an initial vector $\underline{ v}(0) \in {\R}^{n}$, the FLI is provided by
the relation (\cite{FroGuLe00}, \cite{GuLeFro01})
$$
FLI(\underline{ X}(0),  \underline{ v}(0), T) \equiv \sup _{0 < t\leq T}
 \log || \underline{ v}(t)||\, , \qquad  T\geq 0 \  .
$$
 In recent years other numerical tools 
have been introduced: 
the frequency map analysis 
(\cite{Laskar90}, \cite{Laskar92}, \cite{Laskar93}),  
the sup-map analysis (\cite{Laskar94}, \cite{FroLe96a}) 
and the twist angle (\cite{ConVo97},  \cite{FroLe97}).
The twist angle is particularly suited for two dimensional area-preserving
maps. The frequency map  and the sup-map analyses can be applied to a generic system, but they  are more expensive in terms of computational time 
(\cite{FroLe00}).
We remark that the FLI is defined for
each kind of orbit, while the frequency is not well 
defined for chaotic orbits; moreover the FLI provides quantitative 
information about the measure of chaos.
\subsubsection{The standard map as a model problem}
The definition of the FLI can be extended to  discrete systems as follows.
Given a mapping $M$:  ${\R} {^n} \rightarrow {\R}{^n}$, an initial condition
$\underline{ z}(0) \in {\R} {^n}$, and an initial vector $\underline{ v}(0) \in {\R} {^n}$
 of unitary norm, the 
FLI function  is provided by:
$$
FLI(\underline{ z}(0),\underline{ v}(0),T) =\sup _{0 < n\leq T} \log || \underline{ v}(n) || \, , \qquad  T\geq 0 \, ,
$$
where $\underline{ v}(n)$ is defined through the set of coupled equations:
\begin{eqnarray}
\left \{ \begin{array}{lll}
\underline{ z}(n+1) & = & M(\underline{ z}(n)) \nonumber\\
\underline{ v}(n+1) & = &{\partial M\over\partial \underline{ z}}
(\underline{ z}(n))\ \underline{ v}(n) \, .\nonumber
\end{array} \right.
\label{tt} 
\end{eqnarray}
As a model problem, we consider the two dimensional
 standard map (\cite{Fro70}, \cite{LicLie83}, \cite{Chiri79}):
\begin{eqnarray}
M:\left\{ \begin{array}{lllr}
x(n+1) & = & x(n)+\varepsilon 
\sin (x(n)+y(n)) & (\mbox{mod.}\ 2\pi) \\
y(n+1) & = & y(n)+x(n) & (\mbox{mod.}\ 2\pi) \, .
\end{array} \right. 
\label{standard1}
\end{eqnarray}
In Figure \ref{mapp03}, we display some orbits of the standard map defined by 
(\ref{standard1}) with $\varepsilon = 0.6$.
\begin{figure}[htb]
\centering
\includegraphics[height=8cm,width=10cm]{orbi.eps}
\caption{A set of orbits of the standard map for $\varepsilon = 0.6$.}
\label{mapp03}
\end{figure}
For $\varepsilon = 0.6$, the majority of the orbits are
still invariant tori.  Some resonant curves are displayed around 
elliptic points and a 
 chaotic zone 
is generated by the existence of  the hyperbolic point
at the origin. \par
Figure \ref{fli1} shows the variation of the FLI with respect to  time 
for three different kinds of orbits. 
The upper curve, corresponding to initial conditions  $(x,y)=(10^{-4},0)$
selected in the previous chaotic zone,
 shows  an  exponential increase of the FLI; the upper
value of $20$ is a threshold that we impose in order to avoid floating
overflow.
The intermediate curve corresponds to a regular invariant torus 
with initial conditions $(x,y)=(1,0)$ and the lowest one describes a
resonant trajectory with initial conditions $(x,y)=(0,1)$.

\begin{figure}[thb]
\centering
\includegraphics[height=8cm,width=10cm]{flit.eps}
\caption{Variation of the Fast Lyapunov Indicator w.r.t. time for three orbits 
of the standard map with  $\varepsilon = 0.6$. 
The upper curve corresponds to a chaotic
orbit with initial conditions $x(0)=10^{-4}, y(0)=0$; the middle one describes
a non--resonant orbit with $x(0)=1, y(0)=0$ and the lowest one concerns
a resonant orbit with  $x(0)=0, y(0)=1$. }
\label{fli1}
\end{figure}
In \cite{FroLe00} the FLI was computed for a set of 1000 initial 
conditions, regularly spaced on the $x$-axis in the interval $[0,\pi]$,
while $y(0)$ was always set to zero. 
Figure \ref{fli_scan} shows the so called  {\sl FLI-map}, i.e. the value
of the FLI after $T=1000$ iterations against $x(0)$; here, 
$\underline{ v}(0)$  was always taken in the direction of the $x$-axis.
Many orbits appear to have an  $FLI$ equal to the logarithm 
of the integration time, i.e. $\log T = 3$:
actually, this value corresponds to invariant tori
 (\cite{FroLe00}, \cite{GuLeFro01}). \par 
 Values slightly greater than $\log T$ indicate  very thin
chaotic layers or invariant tori close to very thin chaotic zones.
The orbits with an FLI lower than $\log T$ correspond  to chains of islands
(for an explanation see \cite{GuLeFro01}). 
\begin{figure}[bht]
\centering
\includegraphics[height=7cm,width=12cm]{flimap.eps}
\caption{Variation of the FLI as a function 
of the initial
values $x(0)$ for a set of $1000$ orbits regularly spaced in the interval
$[0,\pi]$. The initial vector  $\underline{ v}(0)$ has components $(1,0)$, i.e. it is
almost perpendicular to the corresponding invariant curves.}
\label{fli_scan}
\end{figure}
\subsection{The Arnold's web pictures for a simple Hamiltonian model}
In this section  we give  a high accurate graphical representation of the
 Arnold 
web, obtained implementing the FLI method.
As an example, we consider the following 
Hamiltonian function:
$$
 H_\varepsilon (I_1,I_2,I_3,\varphi _1,\varphi _2,\varphi _3) = {I_1^2 \over 2}+{I_2^2 \over 2} +I_3 + \varepsilon\left(
{ 1 \over \cos (\varphi _1) +
 \cos (\varphi _2) + \cos (\varphi_3) + 4}\right )\ \ ,
$$
where $I_1,I_2,I_3\in {\bf R}$ and $\varphi_1,\varphi_2,\varphi_3 \in {\bf \T}$ are 
canonically conjugated, 
and $\varepsilon$ is a small parameter. The canonical equations of the 
integrable Hamiltonian $H_0$ are trivially integrated: $I_1,I_2,I_3$ stay 
constant, while the angles $\varphi_1(t)=\varphi_1(0)+I_1 t$, $\varphi_2(t)=\varphi_2(0)+
I_2 t$, $\varphi_3(t)=\varphi_3(0)+t$ rotate with 
constant angular velocity. Therefore, each pair of actions
 $(I_1,I_2)$ characterizes an invariant
torus ${\bf T}^3$, on which motions  
are quasi--periodic with frequencies 
$\omega_1=I_1$, $\omega_2=I_2$, $\omega_3=1$. 
Conversely, for any whatsoever small $\varepsilon$ different from zero,
 $H_\varepsilon$ is not integrable. However,  KAM 
theorem provides the existence of a
large volume of invariant tori, 
embedded in the Arnold web. Our goal is to determine 
numerically the structure of the Arnold web, which can be 
 conveniently represented in the 
two--dimensional plane $(I_1,I_2)$. Indeed, each point on this plane 
corresponds to a unique  frequency of an unperturbed torus. Moreover, 
all resonances $k_1\omega_1+k_2\omega_2+k_3\omega_3=0$ are represented by 
the straight lines $k_1I_1+k_2I_2+k_3=0$; the set of all resonances 
is dense on the plane. 
However, one can expect that irregular orbits surround 
each resonance line, up to a distance decreasing as  
$\sqrt{\varepsilon}/ \vert k\vert^\tau$; we refer to this region as
the {\sl resonant zone}. Consequently, the 
volume of the Arnold web is of the order of $\sqrt \varepsilon$. 
We shortly review the qualitative behaviour of the motions with
initial conditions in the Arnold web.
Within resonant zones, both chaotic and regular motions can be observed.
Stable periodic orbits are surrounded by islands on
which the motion is still quasi-periodic, though the dimension of the torus
is strictly smaller than the number of degrees of freedom.
Unstable periodic orbits are surrounded by chaotic 
zones where  Nekhoroshev's theorem  predicts a  diffusion
with a velocity exponentially small with respect
to $-1/\varepsilon$ in the action space. 
When increasing the 
strength of the perturbation, the regular set shrinks until it almost
completely disappears.
In this case the dynamics is no more controlled
by Nekhoroshev's theorem. To describe it, we resort to the well
known Chirikov (\cite{Chiri79}) overlapping criterium,
 which allows the resonant
chaotic orbits to go from one resonance to the other  
(eventually  giving rise to large scale diffusion). 

\subsubsection{From Nekhoroshev to Chirikov}
Using the FLI method, we proceed to describe
the evolution from a mostly ordered system to  a largely chaotic one,
i.e. from the Nekhoroshev to the Chirikov regime.
To this end, we computed the FLI,
 using a leap-frog symplectic
integrator on a grid of $500\times 500$ 
initial conditions regularly spaced in the action space (without
loss of generality, we set
the initial angles  $\varphi_1 =\varphi_2 = \varphi _3 = 0$).
The  initial choice of the 
tangent vector plays a delicate role.
Indeed, it turns out that resonances which are aligned
with the initial tangent vector $(\dot \varphi _1,\dot \varphi_2)$
are not detected by the FLI method.
In order to overcome this problem 
we have chosen  $(\dot \varphi _1,\dot \varphi_2)$, such that
 $\dot \varphi _1/\dot \varphi_2$  is strongly irrational
(the other components $\dot I _1, \dot I _2$
 play a minor role and we set them equal to unity).
Some results are presented in Figure \ref{arno}.
In each  picture, the initial conditions
of $I_1$  and $I_2$ are associated to the corresponding FLI value by a
 different
 colour. 
The lowest values of the FLI appear in black and they correspond to the 
resonant islands of the Arnold web; the highest values appear in yellow
and they correspond to chaotic motions, arising at the crossing nodes of 
resonant lines or arising near the separatrices.
The FLIs corresponding to KAM tori have approximately  the same value 
and therefore they have the same colour (orange in Figure \ref{arno}
 (top, left),
purple in the other figures).
 Therefore, the resonant lines clearly appear to embed large
zones filled by KAM tori 
(see Figure \ref{arno} (top, left) 
and  the enlargement shown in 
Figure \ref{arno} (top, right)). 
  The choice of the colour varies with the value
of the perturbing parameter and with the integration time. 
 Since the perturbation has 
a full Fourier spectrum, i.e. all harmonics are present 
at order $\varepsilon$, a high number of resonances is already present 
for small $\varepsilon$ (Figure \ref{arno} (top));
we remark that all resonances should
appear just by increasing the integration time. 
On the contrary in Figure \ref{arno} (middle), which refers to $\varepsilon = 0.01 $,
the volume of the invariant tori decreases, though the system is still
in the Nekhoroshev regime.  In these figures, following Nekhoroshev theorem,
the chaotic regions become evident at the crossing of the resonances.
In  Figure \ref{arno} (bottom), which refers to 
$\varepsilon = 0.04 $, the dynamical
regime has completely changed. As outlined before, the majority
of invariant tori has disappeared due to  resonance overlapping,
and a big chaotic connected region  replaced the regularity set.
\par
In conclusion,  using a very simple numerical tool, we described the 
structure of the Arnold web  
and its evolution as a function of
the perturbing parameter. We observed a transition 
from an ordered to a chaotic diffusive system, occurring for
$0.01 < \varepsilon _0 < 0.04$ ($\varepsilon _0$ being defined by (\ref{Nek})).
 
\begin{figure}
%\vskip -4 truecm
%\begin{center}
%\setlength{\hoffset}{-5.cm}
\vskip -3 truecm
\includegraphics[height=18cm,width=18cm]{figure3.ps}
%\vskip -6.5 truecm
%\end{center}
%\vskip -4 truecm
\caption{Evolution of the Arnold web for increasing values
of the perturbing parameter $\varepsilon$.
The lowest values of the FLI appear in black for the regular
resonant islands;  the highest values appear in yellow
for   chaotic motions arising at the crossing nodes of 
resonant lines or arising near the separatrix.
The FLI corresponding to KAM tori have about the same value, 
represented in the picture by the same colour 
(orange: top left, purple the others).
 The  colours are adapted to the choice
of $\varepsilon$  and of the integration time. 
 Left column: a large portion of the  action plane;
top:   $\varepsilon=0.001$, $T=1000$; middle: $\varepsilon=0.01$, $T=1000$;
bottom: $\varepsilon=0.04$, $T=1000$. Right column:
enlargement of the figures on the left, obtained with a larger integration 
time; top: $\varepsilon=0.001$, $T=4000$;   
middle: $\varepsilon=0.01$, $T=2000$; bottom:
$\varepsilon=0.04$, $T=2000$ (from \cite{FroGuLe00}). }
\label{arno}
\end{figure}
\subsubsection{ Detection of the diffusion} 
Following Nekhoroshev's theorem, we expect a very slow diffusion of the actions along the resonances, which is very difficult to  detect numerically.
The idea of such kind of diffusion was introduced by the 
pioneering work of Arnold (\cite{Arno64}),  using an ad hoc model. 
Since then, it is called {\sl Arnold's diffusion}. 
 According to (\ref{Nek}), when $\varepsilon$ is close to $\varepsilon _0$, i.e. to
the transition from the stable regime to the diffusive one,
the actions can become quickly unstable.
Therefore, in order to measure Arnold's diffusion, we  need to start 
close to $\varepsilon _0$. 
In the previous section, we found 
that the transition occurred in the interval $0.01 < \varepsilon _0 < 0.04$
(this interval was confirmed and refined to
 $0.0300 < \varepsilon _0 < 0.032$ in \cite{GuLeFro01}).
We used the FLI pictures 
to find chaotic initial conditions on a given 
resonant line as well  as to observe  that diffusion occurs along
this line.  \par
Figure \ref{diffu1} shows some enlargements of the FLI pictures in the action
 space,
around $I_1 = 0.3$, $I_2 = 0.14$ and for different values of $\varepsilon $.
The region between the two yellow lines corresponds to the resonance
associated to $I _1 - 2I_2 = 0$, while the two yellow lines correspond
to its hyperbolic border where diffusion is  confined. The graphical inspection provides the possibility of choosing initial data in the hyperbolic border.
We can follow the evolution of the solution computing the double section 
 of the solution:
  $ \sigma = |\varphi _1| +    
|\varphi _2| \leq 0.05$, $\varphi _3 = 0$. 
This strategy allows to minimize all projection effects and fast quasi-periodic motions; henceforth, we can observe  
a very slow drift along the border of the resonance.\par
We considered  $100$ initial conditions 
 in the  interval 
$0.303 \leq I_1 \leq 0.304$, $0.143 \leq I_2 \leq 0.144$,
having FLI values larger than
$1.2\cdot \log(T)$; 
such initial data  generate chaotic orbits at the border of the resonance and
they are chosen far from the more stable crossing with other resonances.
Let us remark that the points in the double section 
 will appear on both sides of the resonance (in fact  the two yellow lines are connected by a hyperbolic region in the
six dimensional phase space). \par
Figure \ref{diffu1} (top, left) provides the successive intersections with $\sigma \leq 0.05$, $\varphi _3 = 0$,
 up to a time $t = 10^7$; Figure \ref{diffu1} (top, right)
 extends to 
$t=10^8$. 
Even if  the perturbing parameter is one order of magnitude lower than the
 threshold $\varepsilon _0$, 
the  diffusion phenomenon along the resonant line clearly appears.
When decreasing $\varepsilon$ to the minimum
value   $\varepsilon = 0.001$, we detected  diffusion along the resonance
with smaller and smaller speed.
For $\varepsilon=0.007$ (Figure \ref{diffu1} (middle, left)
 for $t=10^6$,  Figure \ref{diffu1} (middle, right) for $t=2.4\cdot 10^7$),  
one obtains similar results, though the speed of diffusion is larger. 
For $\varepsilon = 0.02$  (Figure \ref{diffu1} (bottom, left)
 for $t=1.6\cdot 10^4$, 
 Figure \ref{diffu1} (bottom, right) for $t=5\cdot 10^5$),  diffusion along the resonance is still evident, although it extends a little also on the direction transversal to the resonance. This phenomenon is due to the fact that we are approaching the transition value and  higher order resonances intersecting the main one become evident.
When approaching the critical value, we expect that a large chaotic
 region replaces the zone of invariant tori, giving rise to fast diffusion.
\subsubsection{Measuring  the diffusion coefficient}
Though it does not exist yet an exhaustive analytic model,
we tried to measure a diffusion coefficient as if the phenomenon was Brownian
like.  The numerical experiments are constrained by computational limitations,
since we have to find a compromise between the number of initial conditions
and the length of the integration time. 
Anyway, we could confirm the Brownian hypothesis, i.e. 
we observed an averaged linear increase with time, with slope $D$,
 of the mean squared distance from the initial conditions.
More precisely 
we computed the quantity:
\begin{equation}
S(t) = {1\over N} \sum _{j=1}^N [(I_{2,j}(t)+2I_{1,j}(t))-(I_{2,j}(0)+2I_{1,j}(0))]^2 \, ,
\label{mean}
\end{equation}
 where $I_{1,j}(0)$ and   $I_{2,j}(0)$,
$j=1,...,N$, are 
the initial conditions of a set of $N$  orbits and  
 $I_{1,j}(t)$ and   $I_{2,j}(t)$ are the corresponding values at time $t$.
The term in the squared brackets of (\ref{mean}) represents the  square of the distances
of the actions from the initial values, projected
on the resonant line $I_1 = 2I_2$.
The estimates of $D$ versus   $1/\varepsilon$ are provided in 
Figure \ref{depsilon} on a logarithmic scale. 
Clearly, we can exclude a linear regression,
 which would correspond
to a power law $D(\varepsilon) = C \, (1/\varepsilon)^m$. 
For some  different sets of data,
we performed local regression, finding different slopes.
In particular, the results are the following:
the set containing the values of $D$ for $1/\varepsilon \leq 55 $ has
slope $m = -4.5 $; the set for  $62 \leq 1/\varepsilon \leq 250$  has  
$m = -6.9 $; the set with $1/\varepsilon \geq 330 $  has  $m = -8.8$.
 The changes of the slopes
agree with the expected exponential decrease of $D$.
An exponential
fit of the form $D(\varepsilon) = C'\exp(-\kappa/\varepsilon)^\alpha$ 
(for some parameters  $\kappa$ and $\alpha$) requires
a larger $\varepsilon$--range, which 
is outside our computational possibilities. To give en example,
 the value of $D$ for $1/\varepsilon = 1000$  
(see Figure \ref{depsilon}) required 4 months of CPU time on a 
fast workstation (Compaq AlphaStation XP1000 6/667mhz).

\begin{figure}
\vskip -4 truecm
\includegraphics[height=22cm,width=18cm]{diffu6.ps}
\vskip -1truecm
\caption{ Diffusion along the resonant line $I_1 = 2I_2$ 
for  $\varepsilon = 0.003$ (top), $\varepsilon =0.007$ (middle),
$\varepsilon = 0.02$  (bottom) of a set of $100$ initial conditions  in
the hyperbolic border of the resonance within  the interval 
$0.303 \leq I_1 \leq 0.304$ and $0.143 \leq I_2 \leq 0.144$. 
Black points correspond to the intersections of the orbits on the double 
section 
$|\varphi _1|+|\varphi_2| \leq 0.05$, $\varphi_3 = 0$.
The integration times are respectively: $10^7$ (top, left),
 $10^8$ (top, right), $10^6$ (middle, left), $2.4\cdot 10^7$ (middle, right),
 $1.6\cdot 10^4$ (bottom, left), $5\cdot 10^5$ (bottom, right). The colour
 ranges from black to yellow (from \cite{LeGuFro02}).}
\label{diffu1}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[height=12cm,width=14cm]{slope_new.eps}
\end{center}
\caption{ Computation of the diffusion coefficient as a function of 
$1/\varepsilon$.
The change of the slope of the three power--law fits
agrees with the expected exponential decrease of D
(from \cite{LeGuFro02}).}
\label{depsilon}
\end{figure}


\vglue2cm 

\bf Acknowledgments: \rm One of the authors (A.C.) is deeply indebted to L. Chierchia for 
useful suggestions and discussions. 


\vfill\eject 


\noindent

\begin{thebibliography}{10}
\bibitem{Ar1}V.I. Arnold.  Proof of a Theorem by A.N. Kolmogorov on
the invariance of quasi-periodic motions under small perturbations of the
Hamiltonian. {\sl Russ. Math. Surveys}, {\bf 18}, 9 (1963)

\bibitem{Ar2}V. I. Arnold. { Small denominators and problems of stability of motion 
in classical and Celestial Mechanics}. {\sl Russian Math. Survey},
 {\bf 18}, 85  (1963)

\bibitem{Arno64}
V.I. Arnold.
\newblock Instability of dynamical systems with several degrees of freedom.
\newblock {\em Sov. Math. Dokl.}, {\bf 6}, 581 (1964)
\bibitem{Arbook} V. I. Arnold. {\sl Mathematical methods of classical 
mechanics.}, Springer--Verlag, New York 
(1978) (Russian original, Moscow, 1974) 

\bibitem{Arenc} V. I. Arnold (editor). \sl Encyclopaedia of Mathematical
Sciences. \rm  Dynamical Systems III, Springer-Verlag {\bf 3} (1988)



\bibitem{BeFaGu98}
G.~Benettin, F.~Fass\'o,  M.~Guzzo.
\newblock Nekhoroshev-stability of $L_4$ and $L_5$ in the spatial restricted
  three-body problem.
\newblock {\em Regular and Chaotic Dynamics}, {\bf 3}, 56 (1988)
\bibitem{BeneGa86}
G.~Benettin, G. Gallavotti.
\newblock Stability of motions near resonances in quasi integrable Hamiltonian
systems.
\newblock {\em Journ. Stat. Phys.}, {\bf 44}, 293 (1986)

\bibitem{Czamp}  A. Celletti. Analysis of resonances in the spin--orbit 
problem in Celestial Mechanics: The synchronous resonance (Part I). 
{\sl J. of Appl. Math. and Phys. (ZAMP)}, {\bf 41}, 174 (1990)

\bibitem{CCJMP}  A. Celletti, L. Chierchia. Rigorous estimates for a 
computer-assisted KAM theory. {\sl J. Math. Phys.}, {\bf 28}, 2078 (1987) 

\bibitem{CC0}  A. Celletti, L. Chierchia.  A Constructive Theory
of Lagrangian Tori and Computer--assisted Applications. {\sl Dynamics
Reported} (C.K.R.T. Jones, U. Kirchgraber, H.O. Walther Managing Editors),  
{\bf 4} (New Series), Springer--Verlag, 60 (1995)

\bibitem{CC1}  A. Celletti, L. Chierchia.  { On  the stability of
realistic three-body problems}. {\sl Commun. Math. Phys.},  {\bf 186},
413 (1997)

\bibitem{CC}  A. Celletti, L. Chierchia.  {\it KAM Stability and Celestial Mechanics}. Preprint (2003) 
\bibitem{Chiri79}
B.V. Chirikov.
\newblock A universal instability of many dimensional oscillator system.
\newblock {\em Phys. Rep.}, {\bf 52}, 263 (1979)
\bibitem{ConVo97}
G.~Contopoulos,  N.~Voglis.
\newblock A fast method for distinguishing between order and chaotic orbits.
\newblock {\em Astron. Astrophys.}, {\bf 317}, 73 (1997)
\bibitem{D} C. Delaunay. {\sl Th\'eorie du Mouvement de la Lune}.
M\'emoires de l'Acad\'emie des Sciences {\bf 1}, Tome XXVIII, Paris
(1860)

\bibitem{EW} J.-P. Eckmann, P. Wittwer. {\sl Computer methods and
Borel sommability applied to Feigenbaum's equation}. Springer
Lecture Notes in Physics, {\bf 227} (1985)

\bibitem{FL} C. Falcolini, R. de la Llave.  A rigorous partial
justification of Greene's criterion. {\sl J. Stat. Phys.}, {\bf 67}, n.3/4, 609
(1992)

\bibitem{Fro70}
C.~Froeschl\'e.
\newblock A numerical study of the stochasticity of dynamical systems with two
  degrees of freedom.
\newblock {\em Astron. Astrophys.}, {\bf 9}, 15 (1970)
\bibitem{Fro84}
C.~Froeschl\'e.
\newblock The Lyapunov characteristic exponents and applications.
\newblock {\em Journal de M\'ec. th\'eor. et appl.}, Numero sp\'ecial, 101
  (1984)
\bibitem{FroGuLe00}
C.~Froeschl\'e, M.~Guzzo, E. Lega. 
\newblock Graphical evolution of the Arnold's web,  from order to chaos.
\newblock {\em Science}, {\bf 289} N.5487, 2108 (2000)
\bibitem{FroLe96a}
C.~Froeschl\'e,  E.~Lega.
\newblock On the measure of the structure around the last KAM torus before and
  after its break-up.
\newblock {\em Celest. Mech. and Dynamical Astron.}, {\bf 64}, 21 (1996)
\bibitem{FroLe97}
C.~Froeschl\'e,  E.~Lega.
\newblock Twist angles,  a fast method for distinguishing islands, tori and weak
  chaotic orbits. Comparison with other methods of analysis.
\newblock {\em Astron. Astrophys.}, {\bf 334}, 355 (1998)

\bibitem{FroLe00}
C.~Froeschl\'e,  E.~Lega.
\newblock On the structure of symplectic mappings. The fast Lyapunov 
indicator,   a very sensitive tool.
\newblock {\em Celest. Mech. and Dynamical Astronomy}, {\bf 78}, 167 (2000)
\bibitem{FroLeGo96}
C.~Froeschl\'e, E.~Lega,  R.~Gonczi.
\newblock Fast Lyapunov indicators. Application to asteroidal motion.
\newblock {\em Celest. Mech. and Dynam. Astron.}, {\bf 67}, 41  (1997)

\bibitem{G} G. Gallavotti. {\sl The Elements of Mechanics}. 
Springer--Verlag, New York (1983)  

\bibitem{Gr}  J.M. Greene.  A method for determining a stochastic transition.
{\sl J. of Math. Phys}, {\bf 20}, 1183 (1979)


\bibitem{GuLeFro01}
M.~Guzzo, E.~Lega,  C.~Froeschl\'e.
\newblock On the numerical detection of the stability of chaotic motions in
  quasi--integrable systems.
\newblock {\em Physica D}, {\bf 163}, 1  (2002)
\bibitem{GuMo97}
M.~Guzzo,  A.~Morbidelli.
\newblock Construction of a Nekhoroshev like result for the asteroid belt
  dynamical system.
\newblock {\em Cel. Mech. and Dyn. Astron.}, {\bf 66}, 255 (1997)
\bibitem{He} M. H\'enon. { Explorationes num\'erique du probl\`eme
restreint IV: Masses egales, orbites non periodique}. {\sl Bullettin
Astronomique}, {\bf 3}, no. 1, fasc. 2, 49 (1966)

\bibitem{HenHeil64}
M.~H\'enon,  C.~Heiles.
\newblock The applicability of the third integral of motion: Some numerical
  experiments.
\newblock {\em Astron. J.}, {\bf 1}:73 (1964)

\bibitem{Wi1} H. Koch, A. Schenkel, P. Wittwer. Computer-Assisted
Proofs in Analysis and Programming in Logic: A Case Study. 
{\sl SIAM Review}, {\bf 38}, n. 4, 565 (1996)
 

\bibitem{K} A.N. Kolmogorov.   On the conservation of conditionally
periodic motions under small perturbation of the Hamiltonian.
{\sl Dokl. Akad. Nauk. SSR}, {\bf 98}, 469 (1954)


\bibitem{La}  O.E. Lanford III.  Computer assisted proofs in analysis.
{\sl Physics A},  {\bf 124}, 465 (1984)

\bibitem{Laskar90}
J.~Laskar.
\newblock The chaotic motion of the solar system. A numerical estimate of the
  size of the chaotic zones.
\newblock {\em Icarus}, {\bf 88}, 266 (1990)
\bibitem{Laskar93}
J.~Laskar.
\newblock Frequency analysis for multi-dimensional systems. Global dynamics and
  diffusion.
\newblock {\em Physica D}, {\bf 67}, 257 (1993)
\bibitem{Laskar94}
J.~Laskar.
\newblock Large scale chaos in the solar system.
\newblock {\em Astron. Astrophys.}, {\bf 287}, L9  (1994)
\bibitem{Laskar92}
J.~Laskar, C.~Froeschl\'e,  A.~Celletti.
\newblock The measure of chaos by the numerical analysis of the fundamental
  frequencies. Application to the standard mapping.
\newblock {\em Physica D}, {\bf 56}, 253 (1992)
\bibitem{LeFro96}
E.~Lega,  C.~Froeschl\'e.
\newblock {\em Fast Lyapunov Indicators. Comparison with other chaos
  indicators. Application to two and four dimensional maps.}
\newblock In The Dynamical Behaviour of Our Planetary System. Kluwer academ.
  publ., J. Henrard and R. Dvorak eds. (1997)
\bibitem{LeGuFro02}
E.~Lega, M.~Guzzo,  C.~Froeschl\'e.
\newblock Detection of Arnold diffusion in Hamiltonian systems.
\newblock {\em Physica D}, {accepted}  (2003)
\bibitem{LicLie83}
A.J. Lichtenberg,  M.A. Lieberman.
\newblock {\em Regular and Stochastic motion}.
\newblock Springer, Berlin, Heidelberg, New York (1983)
\bibitem{Lochak92}
P.~Lochak.
\newblock On the conservation of conditionally periodic motions under small
  perturbation of the Hamiltonian.
\newblock {\em Russ. Math. Surv.}, {\bf 47}, 57 (1992)
\bibitem{MK} R.S. MacKay. Greene's residue criterion. 
{\sl Nonlinearity}, {\bf 5}, 161 (1992)

\bibitem{MonthHerrera79}
M.~Month, J.C.Herrera.
\newblock {\em Nonlinear dynamics and the beam-beam interaction.}
\newblock American Institute of Physics, Month and Errera eds. (1979)

\bibitem{M} J. Moser.  On invariant curves of area-preserving
mappings of an annulus.
{\sl Nach. Akad. Wiss. G\"ottingen, Math. Phys. Kl. II}, {\bf 1}, 1 (1962)
\bibitem{Nek77}
N.N. Nekhoroshev.
\newblock Exponential estimates of the stability time of near-integrable
  Hamiltonian systems.
\newblock {\em Russ. Math. Surveys}, {\bf 32}, 1 (1977)

\bibitem{nesMello97}
D.~Nesvorny, S.~Ferraz-Mello.
\newblock On the asteroidal population of the first-order Jovian resonances.
\newblock {\em Icarus}, {\bf 130}, 247 (1997)
\bibitem{PapLas98}
Y.~Papaphilippou, J.~Laskar.
\newblock Global dynamics of triaxial galactic models through frequency map
  analysis.
\newblock {\em Astronomy and Astrophysics}, {\bf 329}, 451 (1998)
\bibitem{Po} H. Poincar\'e. 
{\sl Les m\'ethodes nouvelles de la m\'echanique c\'eleste}. 
Gauthier Villars, Paris (1899) 


\bibitem{Poschel93}
J.~P\"oschel.
\newblock Nekhoroshev's estimates for quasi-convex Hamiltonian systems.
\newblock {\em Math. Z.}, {\bf 213}, 187 (1993)
\bibitem{SM}  C.L. Siegel, J.K. Moser. \sl  Lectures on Celestial Mechanics. \rm
Springer--Verlag Berlin Heidelberg, New York (1971)


\bibitem{Sz} V. Szebehely. {\sl Theory of orbits}. Academic Press, New
York and London (1967)

\end{thebibliography}

\end{document}





