INSTRUCTIONS

The figures of this preprint can be obtained as postscript files with
a request, via computer mail, to one of the authors.


 The text between the lines BODY and ENDBODY is made  of


   816 lines   and       33248 bytes



 In the following table this count is broken down by ASCII  code;
 immediately following the code is the corresponding  character.

     20270  lowercase letters
       747  uppercase letters
       703  digits
       114  ASCII characters   9
       816  ASCII characters  10

      4667  ASCII characters  32
       662  ASCII characters  36 $
        85  ASCII characters  38 &
        73  ASCII characters  39 '
       260  ASCII characters  40 (
       261  ASCII characters  41 )
        20  ASCII characters  42 *
        57  ASCII characters  43 +
       350  ASCII characters  44 ,
       107  ASCII characters  45 -
       254  ASCII characters  46 .
        48  ASCII characters  47 /
        54  ASCII characters  58 :
        17  ASCII characters  59 ;
         6  ASCII characters  60 <
       133  ASCII characters  61 =
         2  ASCII characters  62 >
         3  ASCII characters  63 ?
         4  ASCII characters  64 @
        61  ASCII characters  91 [
      1173  ASCII characters  92 \
        59  ASCII characters  93 ]
        68  ASCII characters  94 ^
       297  ASCII characters  95 _
        31  ASCII characters  96 `
       908  ASCII characters 123 {
        19  ASCII characters 124 |
       905  ASCII characters 125 }
        14  ASCII characters 126 ~

 BODY
\documentstyle[12pt,amssymbols]{article}
\begin{document}
\begin{titlepage}
\addtolength{\baselineskip}{.6\baselineskip}
\noindent\hspace*{.35cm}{\LARGE\bf C}amerino
\hfill {\normalsize CARR REPORTS IN}\\
l' {\LARGE\bf A}quila \hfill {\normalsize MATHEMATICAL PHYSICS}\\
\hspace*{.35cm}{\LARGE\bf R}oma La Sapienza \hfill\\
\hspace*{.35cm}{\LARGE\bf R}oma Tor Vergata \hfill\\
\vfill
    \begin{center}{\large\bf
        Natural Boundaries for Area-Preserving Twist Maps
    }\end{center}
    \vspace{2cm}
    \begin{center}{\sc
A. Berretti, A. Celletti, L. Chierchia, C. Falcolini
    }\end{center}
\vfill
{\normalsize July 1991 \hfill n.10/91}
\end{titlepage}
\citation{CC}
\citation{BC}
\citation{MK1}
\citation{M}
\citation{MK2}
\newlabel{ymap}{{1}{2}}
\newlabel{xmap}{{2}{2}}
\newlabel{dioph}{{3}{2}}
\newlabel{conj}{{6}{3}}
\newlabel{series}{{7}{4}}
\newlabel{cmap}{{8}{4}}
\newlabel{annulus}{{9}{5}}
\newlabel{funct}{{10}{5}}
\citation{BC}
\newlabel{f12}{{11}{6}}
\newlabel{f15}{{12}{6}}
\newlabel{f135}{{13}{6}}
\newlabel{f123}{{14}{6}}
\citation{ED}
\citation{GP}
\citation{Gr}
\citation{Gr}
\citation{Gr}
\newlabel{formalseries}{{15}{11}}
\citation{Baker}
\bibcite{CC}{1}
\newlabel{Frec1}{{16}{13}}
\newlabel{Frec2}{{17}{13}}
\newlabel{Frec3}{{18}{13}}
\bibcite{BC}{2}
\bibcite{MK1}{3}
\bibcite{M}{4}
\bibcite{MK2}{5}
\bibcite{ED}{6}
\bibcite{GP}{7}
\bibcite{Gr}{8}
\bibcite{Baker}{9}

\newcommand{\eps}{\varepsilon}
\newcommand{\om}{\omega}
\newcommand{\th}{\theta}
\newcommand{\ol}{\overline}
    \begin{center}{\Large\sc Natural Boundaries for Area-Preserving}\end{center}
    \begin{center}{\Large\sc Twist Maps}\end{center}
    \vfill
    \begin{center}{\large
        A. Berretti{\small$^{\dag}$},
        A. Celletti{\small$^{\ddag}$},
        L. Chierchia{\small$^{\dag}$},
        C. Falcolini{\small$^{\dag}$}
    }\end{center}
    \vfill
    \begin{center}{\large\bf Abstract}\end{center}

    \noindent We consider KAM invariant curves for  generalizations of
    the standard map of the form $(x',y')=(x+y',y+\eps f(x))$,
    where $f(x)$ is an odd trigonometric polynomial. We study
    numerically their analytic properties by Pad\'e approximant
    method applied to the
    function which conjugates the
    dynamics to a rotation $\theta \mapsto \th + \om$. In the complex
    $\eps$ plane, natural
    boundaries of different shapes are found. In the complex $\th$
    plane the analyticity region appears to be a strip bounded by a
    natural boundary, whose width tends linearly to $0$ as $\eps$
    tends to the critical value.

    \vspace{0.5cm}

{\bf Keywords: } Conservative Dynamical Systems, KAM Theory, Natural
Boundaries, Pad\'{e} Approximants
\vfill
\noindent\setlength{\unitlength}{1mm}
\begin{picture}(100,2)
\put (0,1){\line(1,0){100}}
\end{picture}

\noindent$^{\dag}$Dipartimento di Matematica, II Universit\`a degli
Studi di  Roma (Tor Vergata), Via del Fontanile di Carcaricola, 00133
Roma (Italy).\\
$^{\ddag}$Dipartimento di Matematica Pura e Applicata, Universit\`a
di L'Aquila, 67100 Coppito--L'Aquila (Italy).\\
E-mail: {\tt berretti@roma2.infn.it}, {\tt celletti@vxscaq.infn.it},\\
{\tt chierchia@irmtvm51.bitnet}, {\tt falcolini@irmtvm51.bitnet}.
\pagebreak
\renewcommand{\baselinestretch}{1.5}

\section{Introduction}

The purpose of this paper is to continue and extend the analysis,
 started in \cite{CC} and \cite{BC}, of the complex analytic properties
of invariant curves for area preserving twist diffeomorphims
(for a review, see e.g. \cite{MK1}, \cite{M}).

We consider the following class of area-preserving twist maps
$F_{\eps}: (x,y) \mapsto (x',y')$ of the
cylinder ${\cal C}\equiv \Bbb T \times \Bbb R$ into itself:

\begin{eqnarray}
    y' & = & y + \eps f(x) \label{ymap} \\
    x' & = & x + y'  \pmod{2 \pi}\ , \label{xmap}
\end{eqnarray}

\noindent where $f(x)$ is an odd trigonometric polynomial (therefore with
vanishing mean value) and $\eps$ a real parameter. Such maps
generalize the so-called ``standard map'', where $f(x) = \sin x$.

Given a number $\om$ such that $\om / 2 \pi$ is diophantine, namely
\begin{equation}
\exists \gamma,\tau \geq 1, \mbox{ such that } \forall n \in {\Bbb Z}
\setminus \{0\},
m \in {\Bbb Z}: \left|\frac{\om}{2 \pi}n + m\right| \geq
\frac{1}{\gamma |n|^{\tau}}\ , \label{dioph}
\end{equation}
KAM theory tells that, for $|\eps|$ small enough, there exists
an invariant curve $\Gamma = \Gamma_{\om}(\eps) \subset \cal C$,
homotopically non-trivial (i.e. which winds around the cylinder), invariant
under $F_{\eps}$ ($F_{\eps}(\Gamma_{\om}) = \Gamma_{\om}$) and with
rotation number $\om$. By a theorem of Birkhoff, such curves are graphs of
Lipschitz continuous functions of $x$.

Such ``KAM curves'' have very strong regularity properties. In particular,
for $|\eps|$ small enough, they are analytically conjugated to a rotation by
$\om$; moreover, they are smooth deformations, as $\eps$ grows, of the
trivial invariant curves ($\Gamma_\omega(0)=\{(x,y)\in{\cal C}:y = \om\}$)
obtained at $\eps = 0$.

On the other hand, if $|\eps|$ is large enough, it can be shown (see \cite{MK2})
that no such curves exist at all.
>From now on, by ``invariant curve'' we shall mean a closed, homotopically
non-trivial, invariant curve.

For each $\om \in (0,2\pi)$ we can define two ``thresholds'' in the following
way. Let:
\begin{eqnarray*}
{\cal E}(\om)  & = & \{\eps \geq 0 | \exists
        \Gamma_{\om}({\eps}), \mbox{ jointly analytic in }
        \Bbb T \times [0,{\eps})\}, \\
{\cal E}'(\om) & = & \{\eps \geq 0 | \exists
        \Gamma_{\om}({\eps}), \mbox{ jointly continuous in }
        \Bbb T \times [0,{\eps})\},
\end{eqnarray*}

and:
\begin{eqnarray*}
\eps_{c}(\om)  & = & \sup {\cal E}_\om(\eps), \\
\eps_{c}'(\om) & = & \sup {\cal E}_\om'(\eps).
\end{eqnarray*}

It is believed that, for the standard map, $\eps_{c}(\om) = \eps_{c}'(\om)$,
but no proof of this fact exists. The number $\eps_{c}(\om)$ is usually
called the {\em (analytic) breakdown threshold} for the invariant curve
$\Gamma_{\om}$. As a function of $\om$, $\eps_{c}(\om)$ is quite irregular,
being in general 0 for each $\om$ rational and non-zero for $\om$ diophantine.

The breakdown mechanism and, in general, the properties of invariant curves
near the critical threshold $\eps_{c}(\om)$ are far to be understood.

Here we shall investigate numerically the regularity properties of
the invariant curves and in particular their analytic structure.

First of all, we note that the dynamics may be described in terms
of the variable $x \in \Bbb T$ only; in fact, $y' = x' - x$ and therefore:
\begin{equation}
x_{n + 1} - 2 x_{n} + x_{n - 1} = \eps f(x_{n}),
\end{equation}
where $x_{n} = F_{\eps}^{n}(x_{0},y_{0})$.
Invariant curves may then be found by looking for a change of variables
in which the dynamics reduces to a simple rotation by $\om$:
\begin{equation}
x = \th + u(\th),
\end{equation}
and:
\[
\th_{n} = \th_{n - 1} + \om = n \om + \th_{0},
\]
where of course it must be $1 + u'(\th) > 0$. It is simple to check that
the invariant curve with rotation number $\om$ is given, in parametric form, by:
\begin{eqnarray*}
x & = & \th + u(\th) \\
y & = & \om +u(\th) - u(\th - \om)
\end{eqnarray*}
and that the invariant curve shares the same regularity properties of
$u(\th)$, by the implicit function theorem and the above mentioned theorem of
Birkhoff.

The function $u(\th)$ (from now on called briefly ``conjugating function'')
satisfies the following equation:
\begin{equation}
D_{\om}^{2}u \equiv u(\th + \om) - 2 u(\th) + u(\th - \om) =
        \eps f(\th + u(\th)). \label{conj}
\end{equation}
This equation may be studied perturbatively by expanding $u(\th)$ in
powers of $\eps$ and, further, the Taylor coefficients of $u(\th)$ as
Fourier series in $\th$:
\begin{equation}
u(\th;\eps) = \sum_{n = 1}^{\infty} u_n(\th)\eps^n =
\sum_{n = 1}^{\infty} \sum_{k \in {\Bbb Z}}
        \hat{u}_{n,k} e^{i k \th} \eps^{n} .    \label{series}
\end{equation}
Actually, since in our case $f(x)$ is always a trigonometric polynomial
containing only sines, $u(\th)$ is odd and $\hat{u}_{n,k} \in i {\Bbb R}$,
$\hat{u}_{n,-k} = - \hat{u}_{n,k}$, as it is very easy to check.
Moreover, each Taylor coefficient will be a trigonometric polynomial
of order increasing with $n$.

It is natural to study the analytic properties of $u(\th)$ by investigating
the series (\ref{series}). In particular, we are interested {\em both}
in the analyticity in $\eps$ at each fixed $\th$, {\em and} the analyticity
in $\th$ at fixed $\eps$. In particular, a natural question which arises is
to find
the relation between the radius of convergence in the complex $\eps$ plane
of (\ref{series}) at fixed $\th$ and the KAM break-down threshold defined
above.
To be precise, let:
\[
\rho(\om) = \inf_{\th \in \Bbb T} \left( \limsup_{n \rightarrow \infty}
        |u_{n}(\th)|^{1/n} \right)^{-1}
\]
be the radius of convergence, uniform in $\th$, of the series (\ref{series}).
Clearly $\rho(\om) \leq \eps_{c}(\om)$, since within the radius of
convergence of the series (\ref{series}) there exists an analytic
solution to the equation (\ref{conj}).
It is interesting to study the exact relation between this two
``thresholds'', as well as to understand the mechanism by which the
series (\ref{series}) ceases to converge as $|\eps|$ grows and
therefore $u(\th;\eps)$ looses analyticity in $\eps$.

A different but related question is the behavior of the same series
expansion on the complex $\th$ plane
at {\em fixed $\eps$ within} the radius of convergence.
It is more convenient to pass to
complexified dynamical variables as follows:
\begin{eqnarray*}
z     & = & e^{i x}, \\
\zeta & = & e^{i \th}, \\
P(z)  & = & i f(x).
\end{eqnarray*}
With this notations, the maps we consider can be written as:
\begin{equation}
\frac{z_{n+1} z_{n-1}}{z_{n}^{2}} = e^{\eps P(z_{n})} \label{cmap}
\end{equation}
where $P(z)$ is a rational function of $z$ with a pole at the origin,
such that the coefficients $\{c_{n}\}$ of its Laurent series at $z=0$
are real and satisfy $c_{-k} = - c_{k}$, $c_{0} = 0$. In particular,
the unit circle in the complex $z$ plane is invariant under this
complex dynamics.

Next, let $\phi(\zeta) = e^{i (\th + u(\th))}$ and
$\gamma = e^{i \om}$; then $\phi$ conjugates multiplication by
$\gamma$ (i.e. rotation by $\om$) on the $\zeta$ plane to the
dynamics of our complexified map on the $z$ plane. In particular,
circles around the origin of radius close enough to 1 on the
$\zeta$ plane are mapped conformally to analytic, closed invariant
curves which wind around the unit circle, as a result of the analytic
KAM theory.

The series (\ref{series}) can be written as a double Taylor-Laurent
series:
\[
\sum_{n=1}^{\infty} u_{n}(\th) \eps^{n},
\]
with:
\[
u_{n}(\th) = \sum_{k \in \Bbb Z} {\hat u_{n,k}} \zeta^{k},
\]
and, for each fixed $\eps < \rho({\om})$, the Laurent series in $\zeta$
will converge in an annulus:
\begin{equation}
\sigma^{-1}(\om,\eps) < |\zeta| < \sigma(\om,\eps), \label{annulus}
\end{equation}
the two radii being reciprocal for obvious symmetry reasons.
We investigate the analytic structure of this series as well, and
moreover we study the behavior of $\sigma(\om,\eps)$ as $\eps
\rightarrow \rho(\om)$. Clearly the annulus in the $\zeta$ plane
translates to a strip in the complexified $\th$ plane.

\section{The models}

We consider maps with:
\begin{equation}
f(x) = \sum_{\nu = 1}^{\bar{\nu}} \alpha_{\nu} \sin \nu x, \label{funct}
\end{equation}
where $\alpha_{\nu} \in {\Bbb R}$, $\bar{\nu} \in {\Bbb N}$ and we will take
$\alpha_{1} = 1$.

In particular, we analyzed in detail the following two- or
three-frequencies maps:
\begin{equation}
f_{1}(x) = \sin x + \frac{1}{20} \sin 2 x\ , \label{f12}
\end{equation}
\begin{equation}
f_{2}(x) = \sin x + \frac{1}{50} \sin 5 x\ , \label{f15}
\end{equation}
\begin{equation}
f_{3}(x) = \sin x + \frac{1}{30} \sin 3 x + \frac{1}{50} \sin 5 x\,\label{f135}
\end{equation}
\begin{equation}
f_{4}(x) = \sin x + \frac{1}{20} \sin 2 x + \frac{1}{30} \sin 3 x\.\label{f123}
\end{equation}
Next, we choose a diophantine (\ref{dioph}) rotation number $\om$.
In order to fix the notation let us introduce the {\em continued
fraction} expansion.
For any irrational number $\om \in (0,1)$ let $\{a_k\}_{k\in
{\Bbb N}}$ be the sequence of integer numbers such that
\[
{\frac{\om}{2\pi}} = \frac{1}{a_{1} + \frac{1}{a_{2}+ \frac{1}{a_{3}+\ddots}}};
\]
symbolically we write:
\[
{\frac{\om}{2\pi}} = [a_{1},a_{2},a_{3},\ldots].
\]
A trivial computation shows that the continued fraction expansion of the
golden ratio ${\frac{\om}{2\pi}}\equiv ({\sqrt{5}}-1)/2$ is composed only
by one:
\[
{\frac{\om_1}{2\pi}} = [1,1,1,\ldots]\ = [1^{\infty}].
\]
In analogy to $\om_{1}$ we shall consider the numbers
\[
{\frac{\om_2}{2\pi}} = [2,2,2,\ldots] = [2^{\infty}]
\]
and:
\[
{\frac{\om_3}{2\pi}} = [3,3,3,\ldots] = [3^{\infty}]
\]
(sometimes referred to as {\em silver} and {\em bronze} numbers,
respectively). Beside $\om_{1}$, $\om_{2}$ and $\om_{3}$ we
consider also {\em noble} numbers (i.e. whose continued fraction is
definitely one) of the form:
\[
{\frac{{\tilde\om_{k}}}{2\pi}}= [1,k,1,1,\ldots] = [1,k,1^{\infty}].
\]

\section{Results}

For all the maps considered we found natural
boundaries in the complex $\eps$ plane, as it
was found for the standard map in \cite{BC}.
The natural boundaries appear to be closed curves, symmetric with
respect to the real axis (since the coefficients of the Taylor
series in (\ref{series}) are real) but not circles.

A simple symmetry argument shows that when the Fourier expansion of $f(x)$
contains only odd frequencies, then
the natural boundary must be symmetric with respect to the
imaginary $\eps$ axis (cfr. fig. 2,3): in fact, if $f(x)$ contains only odd
frequencies, $P(z)$ is an odd function and (\ref{cmap}) is invariant
under the transformation $z \mapsto -z$, $\eps \mapsto -\eps$. This
symmetry property implies the observed symmetry in the distribution of the
singularities of $u$.

The intersection of the natural boundary with the positive real axis in
the complex $\eps$ plane coincides with the threshold found by
Greene's method, within numerical errors
(see table 1; a review of Greene's method is presented in App.~A).
It is interesting to note that, since for some maps
the natural boundary appears to elongate in the direction of the real
axis (see table 1 and fig. 2, where it is shown the natural boundary for
{\it e.~g.} the invariant curve with rotation number
$\om_{3}$ of the map with nonlinear term $f_{2}$),
it follows that for those maps the radius of convergence
$\rho(\om)$ of the series (\ref{series}) is {\em strictly less} than
the breakdown threshold $\eps_{c}(\om)$.

We note also that for all the maps of the type we considered (including
those of which we do not show the singularities of Pad\'{e}
approximants for brevity), for all
diophantine rotation numbers of the above mentioned type, a natural
boundary in the complex $\eps$ plane was found. Some natural questions
therefore arise: are there maps with a sufficiently regular nonlinear
term $f$ such that a different behavior occurs? How ``universal'' natural
boundaries are in perturbation expansions for invariant tori of Hamiltonian
systems? What is the relation passes between the shape of the natural
boundaries,
the frequencies in the Fourier expansion of $f$ and the continued fraction
expansion of $\om$?
If indeed we face a universal feature -- at least within the class
of maps considered -- it should be possible to understand it in a
``renormalization group'' framework of the type devised in \cite{ED}.


In the complex $\zeta$ plane, all the maps we considered have natural
boundaries on the annulus of convergence of the Laurent series
(\ref{annulus}).
A similar natural boundary was found, with very different techniques, by
Greene and Percival in \cite{GP} for the standard map and for the
semi-standard map. We show the same boundary directly by showing the
singularities of Pad\'{e} approximants of higher order in the complex
$\zeta$ plane. Moreover, the same phenomena occurs for all the other maps
considered, remarkably always with boundaries of annular shape -- quite
differently than the case of boundaries on the complex $\eps$ plane, where
the shape depends on $f$ and on $\om$.

We also measured the external radius of the annulus
$\sigma(\om,\eps)$ and determined numerically its dependence on $\eps$. In
the case of the semi-standard map ({\it i.~e.} the case $f(x) = e^{i x}$),
since the order $n$ in $\eps$ of the
expansion (\ref{series}) contains only the frequency $k = n$, the expansion can
be considered one in the variable $\eps e^{i \th} = \eps \zeta$; it follows
that $\sigma(\om,\eps)$ tends to $0$ linearly as $\eps$ tends to the radius
of convergence of the series. In the case of the standard map and of the
other more complex maps considered in this paper, this argument of course
does not work, but it is remarkable to note that the width of the
analyticity region in the $\zeta$ plane still appears to tend linearly to
$0$ (and $\sigma(\om,\eps)$ to $1$) as $\eps \rightarrow \eps_{c}(\om)$.

Again, this phenomenology seems to be ``universal'' among the maps
considered here. It is rather easy, though, to write a counter-example
in which {\em no} natural boundary appears for a given invariant curve at a
given $\eps$; in particular, we can exhibit a map of the type
(\ref{ymap}), (\ref{xmap}) which has an invariant curve, say with
rotation number equal to the golden mean, such that its conjugating
function is an {\em entire} function of $\th$ for {\it e.~g.}
$\eps = 1/2$. In fact, let:
\[
u_{\om_{1}}(\th,\frac{1}{2}) = \frac{1}{2} \sin \th
\]
and let $w(x)$ be the inverse function of $\th + 1/2 \sin \th$; then, since for
{\em all} $\om$ diophantine, $\eps$ and $\th$ the equation (\ref{conj})
must be satisfied, a simple calculation shows that if we take:
\[
f(x) = 2 (\cos \pi (\sqrt{5} - 1) - 1) \sin w(x)
\]
then (\ref{conj}) will have $1/2 \sin \th$ as solution for $\om = \om_{1}$
and $\eps = 1/2$.

Besides the examples considered in this paper, we studied a large variety
of mappings of the type given by eq. (\ref{ymap}), (\ref{xmap}),
(\ref{funct}) with rotation numbers of the type mentioned in sect. 2. The
results (not presented here for brevity) show a phenomenology similar to the
cases studied here. We just mention some of them in table 1, where the
radius of convergence and the breakdown threshold as determined by the
Pad\'{e} approximants method and the breakdown threshold as determined with
Greene's method are reported.

Some selected figures showing the shapes of the natural boundary in the
complex $\eps$ plane (fig. 1--4) and in the complex $\zeta$ plane (fig.
5--7) are included.

Finally, in fig. 8 and 9 we show how $\sigma(\om,\eps)$ tends to 1
{\em linearly} as $\eps \rightarrow \rho({\om})$.


\appendix\section{Greene's method}

Up to now, the most reliable numerical method to determine the breakdown of
invariant curves is the one developed by Greene in \cite{Gr}.
The results provided by this method will form the basis of comparison
and interpretation of ours. In this perspective we consider
worthwhile to dedicate a few words to a review of the method.

Greene's starting point is the conjecture that the breakdown of an invariant
curve is related to a transition from stability to instability of the
periodic orbits approaching the invariant curve. More precisely, for any
irrational $\om$ let $\{p_{k}/q_{k}\}_{k \in {\Bbb Z}}$
be the sequence of rational approximants to $\om$ such that
$p_{k}/q_{k} \rightarrow \om$ as $k \rightarrow \infty$.
Let us denote by ${\cal P}(p_{k}/q_{k})$ a periodic orbit with frequency
$p_{k}/q_{k}$. The critical break-down threshold is found as the value of
$\eps$ at which the periodic orbits ${\cal P}(p_{k}/q_{k})$ are {\em
alternatively} stable or unstable. To determine the stability of
${\cal P}(p_{k}/q_{k})$ we compute
the Floquet multipliers of the linearized map. In particular,
the tangent space orbit $(\delta y_{n},\delta x_{n})$ at the point
$(y_{n},x_{n})$ is defined in terms of the initial conditions
$(\delta y_{0},\delta x_{0})$ through a matrix $A$ as:
\[
\left[\begin{array}{c}
        \delta y_{n} \\
        \delta x_{n}
\end{array}\right]      = A     \left[\begin{array}{c}
                                        \delta y_{0} \\
                                        \delta x_{0}
                                \end{array}\right],
\]
where, in our case, the matrix $A$ associated to the periodic orbit
${\cal P}({p_{k}/q_{k}})$ is given by
\[
A = \prod_{i = 1}^{q_{k}}
        \left[\begin{array}{cc}
                1 - \eps \sum_{\nu = 1}^{\bar{\nu}}
                        \alpha_{\nu} \nu \cos \nu x_{i} & 1 \\
                - \eps \sum_{\nu = 1}^{\bar{\nu}}
                        \alpha_{\nu} \nu \cos \nu x_{i} & 1
        \end{array}\right].
\]

The Floquet multipliers of the linearization are the eigenvalues
$\lambda_{1,2}$ of the matrix $A$. Since the map is area-preserving,
such eigenvalues depend only on the trace of $A$.
More precisely, if we set:
\[
R = \frac{1}{4} (2 - \mbox{Trace} (A)),
\]
(called by Greene the {\em residue} of the periodic orbit),
the eigenvalues $\lambda_{1,2}$ of $A$ are given by the expression:
\[
\lambda_{1,2} = 1 - 2 R \pm 2 [R (R - 1)]^{1/2}.
\]

When $0 < R < 1$ the periodic orbit is stable, since the eigenvalues are
complex with modulus $1$. On the other hand the orbit is unstable
when $R < 0$ or $R > 1$.

Therefore once a periodic orbit is found, the calculation of the
residue $R$ determines its stability. The breakdown threshold of an
invariant curve (with rotation number $\om$) is then defined as the
value of $\eps$ at which the periodic orbits ${\cal P}({p_{k}/q_{k}})$
(with frequencies equal to the rational approximants to $\om$)
are alternatively stable and unstable.

Numerically, the most delicate point is certainly the determination
of the periodic orbits. By definition of periodic orbit, to find
${\cal P}({p_{k}/q_{k}})$ one has to find a point $(x_{0},y_{0})$
which is equal to the $q_{k}$-th iterate of the mapping
($x_{0} = x_{q_{k}}$, $y_{0} = y_{q_{k}}$),
with the further constraint that $\sum_{i = 1}^{q_{k}} y_{i} = p_{k}$.
Here $(x_{q_{k}},y_{q_{k}})$ are obtained by the recursive relations:
\[
\left\{
        \begin{array}{ccccl}
                y_{n + 1} & = & y_{n} & + & \eps \sum_{\nu = 1}^{\bar{\nu}}
                        \alpha_{\nu} \sin \nu x_{n}\\
                x_{n + 1} & = & x_{n} & + & y_{n + 1}
        \end{array}
\right.
\]

Notice that since (by definition) $q_{k} \rightarrow \infty$ as
$k \rightarrow \infty$, the determination of ${\cal P}({p_{k}/q_{k}})$
may require a substantial amount of CPU time for large $k$.

To ease this problem, Greene observed (see Appendix $A$ of \cite{Gr}),
using the symmetry of the map, that the initial point of the periodic
orbit belongs to the lines $x_{0} = y_{0}/2$, $y_{0}/2 + \pi$
or $x_{0} = 0$, $\pi$.
In fact, the map can be written as a product of two involutions:
\[
F_{\eps} = I_1 I_2
\]
with $I_{1}^2=I_{2}^2=1$; it is then possible to show \cite{Gr}
that if $(x_{0},y_{0})$ is a fixed point of $I_{1}$ (or $I_{2}$)
and $F_{\eps}^{N}(x_{0},y_{0})$ is also a fixed point for some
$N \in {\Bbb N}$, then the orbit with initial point $(x_{0},y_{0})$
is periodic with period $2 N$.

A simple calculation shows that $I_1$ and $I_2$ are given by:
\[ I_{1}: \left\{
             \begin{array}{ccl}
                x_{n} & = & - x_{n - 1} \\
                y_{n} & = &
                   y_{n - 1} + \eps \sum_{\nu=1}^{\bar{\nu}}
                   \alpha_{\nu} \sin \nu x_{n}
             \end{array}
          \right.
\]
and:
\[ I_{2}: \left\{
             \begin{array}{ccl}
                x_{n + 1} & = & -x_{n} + y_{n} \\
                y_{n}     & = & y_{n}
             \end{array}
          \right. .
\]

The fixed points of $I_{1}$ are at $x = 0$ and $x = \pi$, while those of
$I_{2}$ are at $x = y/2$ or $x = y/2 + \pi$.

\section{Calculation of Pad\'{e} approximants}

Given a formal Taylor series:
\begin{equation}
\sum_{n = 0}^{\infty} a_{n} z^{n}, \label{formalseries}
\end{equation}
its Pad\'{e} approximants are the rational approximants,
traditionally denoted by $[M,N]$, defined in the following way:
\[
[M,N] \equiv \frac{P_{M}(z)}{Q_{N}(z)}
\]
where $P_{M}(z)$, $Q_{N}(z)$ are polynomials of degree $M$, $N$
respectively, such that the Taylor expansion of their quotient agrees up
to order $M + N$ with the series (\ref{formalseries}). $Q_{N}(z)$ is
usually normalized by the condition $Q_{N}(0) = 1$. In the Pad\'{e}
approximant $[M/N]$ there are therefore $M + N + 1$ undeterminate
coefficients as in any polynomial of degree $M + N$.
These coefficients can be formally determined from
first $M + N + 1$ terms of (\ref{formalseries}).
To minimize round off errors, the Pad\'{e} are computed
recursively with the following formulas:
\begin{eqnarray*}
\frac{P_{2j}(x)}{Q_{2j}(x)} & = & [N - j/j] \\
\frac{P_{2j + 1}(x)}{Q_{2j + 1}(x)} & = & [N - j - 1/j],
\end{eqnarray*}
and:
\begin{eqnarray*}
\frac{P_{2j}(x)}{Q_{2j}(x)} & = &
    \frac{(\bar{P}_{2j - 1} P_{2j - 2}(x) - x \bar{P}_{2j - 2} P_{2j - 1}(x))/
        \bar{P}_{2j - 1}}
    {(\bar{P}_{2j - 1} Q_{2j - 2}(x) - x \bar{P}_{2j - 2} Q_{2j - 1}(x))/
        \bar{P}_{2j-1}} \\
\frac{P_{2j + 1}(x)}{Q_{2j + 1}(x)} & = &
    \frac{(\bar{P}_{2j} P_{2j - 1}(x)- \bar{P}_{2j - 1} P_{2j}(x))/
        (\bar{P}_{2j} - \bar{P}_{2j - 1})}
    {(\bar{P}_{2j} Q_{2j - 1}(x)- \bar{P}_{2j - 1} Q_{2j}(x))/
        (\bar{P}_{2j} - \bar{P}_{2j - 1})},
\end{eqnarray*}
where $\bar{P}_{j}$ is the coefficient of the higest power in $P_{j}(x)$,
and the recursion is initialized by setting $P_{0}(x)$ and $P_{1}(x)$
to the $N$ and $N-1$ Taylor polynomial respectively:
\begin{eqnarray*}
P_{0}(x) & = & \sum_{k=0}^{N} a_{k} x^{k}, \mbox{ } Q_{0}(x)= 1, \\
P_{1}(x) & = & \sum_{k=0}^{N-1} a_{k} x^{k}, \mbox{ } Q_{1}(x)= 1.
\end{eqnarray*}
For a complete reference see \cite{Baker}.

To identify a natural boundary we compute ``diagonal'' Pad\'{e}
approximants $[N/N]$ (considered to be the most accurate)
for increasing $N$ and look at how the poles of the
approximants behave: in presence of a natural boundary, they will cluster
on the boundary curve as $N$ grows. For a reliable identification it is
necessary in our case to compute Pad\'{e} approximants up to high orders,
say $[70/70]$ or $[80/80]$.

It is common to have, besides ``genuine'' poles and zeros, also fake
pole-zero pairs which cancel (a so called ``ghost'').
It is possible to distinghish them from a genuine pole and a
genuine zero nearby (a situation very frequent when dealing with
natural boundaries, essential singularities and accumulation
points of  singularities) because ghosts disappear as any parameter
in the series is slightly varied, or the order of the approximants changes.
Moreover, the distance between the pole and the zero in a ghost is always
close to machine precision, and in any case several orders of magnitude
smaller than the distance in a genuine pole-zero pair.
A semi-automatic filtering mechanism can therefore be used.
Finally we remark that since the series we consider have an almost
lacunar nature, as $N$ increases ghosts tend to appear and disappear
with a certain regularity as peaks are crossed.

\section{Calculation of the perturbation expansion}

The coefficients of the perturbative expansion (\ref{series}) are computed
with the following formulas, valid when $f(x)$ is given by (\ref{funct}):
\begin{eqnarray*}
b_{0}^{(\nu)}(\th) & = & e^{i \nu \th} \mbox{\hspace{5mm} for } \nu = 1,
        \ldots, \bar{\nu} \\
u_{n}(\th) & = & D_{\om}^{-2} {\rm Im} \left(\sum_{\nu = 1}^{\bar{\nu}}
\alpha_{\nu} b_{n - 1}^{(\nu)}\right) \mbox{\hspace{5mm} for } n \geq 1, \\
b_{n}^{(\nu)}(\th) & = & \frac{i \nu}{n} \sum_{l = 1}^{n} l u_{l}(\th)
        b_{n - l}^{(\nu)}(\th) \mbox{\hspace{5mm} for } n \geq 1,
\end{eqnarray*}

Next we expand in Fourier series $u$ and the $b$'s; first we note that
since $u(\th)$ is real and odd in $\th$, $\hat{u}_{k}$ is purely imaginary
and odd in $k$, so we set $\hat{v}_{n,k} = i \hat{u}_{n,k}$, with
$\hat{v}_{n,k} \in {\Bbb R}$. We obtain:
\begin{eqnarray}
\hat{b}_{0,k}^{(\nu)} & = & \delta_{\nu,k} \label{Frec1} \\
\hat{v}_{n,k}^{\nu} & = & \frac{1}{2 \hat{D}_{\om,k}^{2}}
        (\sum_{\nu = 1}^{\bar{\nu}} \alpha_{\nu}
        (\hat{b}_{n-1,k}^{(\nu)} - \hat{b}_{n-1,-k}^{(\nu)}) \label{Frec2} \\
\hat{b}_{n,k}^{(\nu)} & = & \frac{\nu}{n} \sum_{l = 1}^{n} l
        \sum_{h = h_{min}}^{h_{max}}
        \hat{v}_{l,h} \hat{b}_{n - l,k - h}^{(\nu)} \label{Frec3}
\end{eqnarray}
with $n \geq 1$; in (\ref{Frec2}) $1 \leq k \leq \bar{\nu} n$
while in (\ref{Frec1}) and (\ref{Frec3})
$- \bar{\nu} n \leq k \leq \bar{\nu} n$; moreover
$h_{min} = - \min (l, (n - l) - k)$,
$h_{max} = \min (l, (n - l) + k)$ and
$\hat{D}_{\om,k}^{2} = 2 (\cos k \om - 1)$.

\section*{Acknowledgements}

All computations have been performed on several VAX/VMS computer systems at
the Dipartimento di Matematica, II Universit\`{a} di Roma (Tor Vergata),
Facolt\`{a} di Scienze, Universit\`{a} dell'Aquila, and Institute of
Scientific Interchange, Torino.
We are grateful to R. De La Llave and G.
Gallavotti for very helpful suggestions and encouragement.
A.~B. and A.~C. whish to express gratitude
to Prof. M. Rasetti, R. Livi and S. Ruffo for their invitations at the
workshops ``Complexity and Evolution'' at the Institute of
Scientific Interchange, Torino, where part of this work was done, and for
many helpful discussions.


\begin{thebibliography}{99}

\bibitem{CC}
Celletti, A., Chierchia, L.: {\it Construction of analytic KAM surfaces and
effective stability bounds}, Commun. in Math. Physics {\bf 118}, 119 (1988)

\bibitem{BC}
Berretti, A., Chierchia, L.: {\it On the complex analytic structure of
the golden invariant curve for the standard map}, Nonlinearity {\bf 3},
39-44 (1990)

\bibitem{MK1}
MacKay, R.~S.: {\it Transition to Chaos for Area-preserving Maps},
Springer Lect. Notes in Physics, vol. 247, p. 390 (1985)

\bibitem{M}
Moser, J.: {\it Recent Developments in the Theory of Hamiltonian
Systems}, SIAM Review {\bf 28}, 459 (1986)

\bibitem{MK2}
Mather, J.~N., Ergodic Th. and Dyn. Sys. {\bf 4}, 301 (1984);
MacKay, R.~S., Percival, I.~C.: {\it Converse KAM: Theory and Practice},
Commun. in Math. Physics {\bf 98}, 469 (1985)

\bibitem{ED}
Escande, D.~F., Doveil F.: {\it Renormalization method for computing the
threshold of large-scale stochastic instability in two degrees of freedom
Hamiltonian systems}, Journal of Stat. Phys. {\bf 26}, 257 (1981)

\bibitem{GP}
Greene J.~M, Percival I.~C., {\it Hamiltonian maps in the complex plane},
Physica {\bf 3D}, 530 (1981)

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


\bibitem{Baker}
Baker G., {\it Essentials of Pad\'{e} Approximants}, Academic Press,
New York 1975

\end{thebibliography}

\newpage

\begin{center}{\bf Table 1}\end{center}
\vspace{2cm}
\begin{center}
\begin{tabular}{c|lll}
                             & $\ol{\rho}$ & $\ol{\eps}_{c}$ &
                                        $\ol{\ol{\eps}}_{c}$ \\ \hline
$f_{1}(x)$, $\om_{1}$         & $0.8$       & $1.2$           &
                                        $1.2$  \\
$f_{2}(x)$, $\om_{3}$         & $0.6$       & $0.7$ &
                                        $0.67$ \\
$f_{3}(x)$, $\om_{1}$         & $0.5$       & $0.5$ &
                                        $0.5$  \\
$f_{4}(x)$, $\om_{1}$         & $0.7$       & $0.9$ &
                                        $0.9$  \\
$f_{1}(x)$, $\tilde{\om_{1}}$ & $0.75$      & $1.1$ &
                                        $1.13$ \\
$f_{2}(x)$, $\tilde{\om_{1}}$ & $0.45$      & $0.6$ &
                                        $0.64$ \\
\end{tabular}
\end{center}
\vspace{2cm}

\noindent
Here $\ol{\rho}$ is an estimate of the radius of convergence of the series
(\ref{series}) based on Pad\'{e} approximants; $\ol{\eps_{c}}$ is an
estimate of the breakdown threshold based on Pad\'{e} approximants;
$\ol{\ol{\eps_{c}}}$ is an estimate of the breakdown threshold based on
Greene's method. The maps $f_{1}$, \ldots, $f_{4}$ have been defined in
sect. 2, eq. (\ref{f12}), \ldots, (\ref{f123}) and the notation for the
rotation numbers is the same as in sect. 2.

\newpage

\begin{center}{\bf Figure Captions}\end{center}
\vspace{2cm}
\begin{description}
\item[Fig. 1] Poles of the Pad\'{e} approximant $[70/70]$ in $\eps$ for
        the map with
        $f(x) = \sin x + 1/20 \sin 2 x$, rotation number $\om_{1}$, $\th = 1.$
\item[Fig. 2] Poles of the Pad\'{e} approximant $[70/70]$ in $\eps$
        for the map with
        $f(x) = \sin x + 1/50 \sin 5 x$, rotation number $\om_{3}$, $\th = 1.$
\item[Fig. 3] Poles of the Pad\'{e} approximant $[70/70]$ in $\eps$
        for the map with
        $f(x) = \sin x + 1/30 \sin 3 x + 1/50 \sin 5 x$,
        rotation number $\om_{1}$, $\th = 1.$
\item[Fig. 4] Poles of the Pad\'{e} approximant $[70/70]$ in $\eps$
        for the map with
        $f(x) = \sin x + 1/20 \sin 2 x + 1/30 \sin 3 x$,
        rotation number $\om_{1}$, $\th = 1.$
\item[Fig. 5] Poles of the Pad\'{e} approximant $[80/80]$ in $\zeta$ for
        the standard map, $\eps = 0.7$.
\item[Fig. 6] Poles of the Pad\'{e} approximant $[80/80]$ in $\zeta$ for
        the standard map, $\eps = 0.8$.
\item[Fig. 7] Poles of the Pad\'{e} approximant $[80/80]$ in $\zeta$ for
        the standard map, $\eps = 0.9$.
\item[Fig. 8] $\sigma(\omega,\eps)$ vs. $\eps$ for the standard map,
        rotation number $\om_{1}$.
\item[Fig. 9] $\sigma(\omega,\eps)$ vs. $\eps$ for the map with
        $f(x) = \sin x + 1/50 \sin 5 x$, rotation number $\om_{1}$.
\end{description}

\end{document}
ENDBODY
