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}