BODY: \newcommand{\Real} {\Bbb R} \newcommand{\Integer} {\Bbb Z} \newcommand {\Natural} {\Bbb N} \newcommand {\Rational} {\Bbb Q} \newcommand{\Complex} {\Bbb C} \newcommand{\Torus} {\Bbb T} \newcommand{\eq}[1]{eq.~(\ref{#1})} %Physica D style \newcommand{\Eq}[1]{Equation~(\ref{#1})} %Physica D style \newcommand{\eqs}[1]{eqs.~(\ref{#1})} %Physica D style \renewcommand{\S}[1]{Sec.~\ref{#1}} %Physica D style \newcommand{\fig}[1]{fig.~(\ref{#1})} %Physica D style \newcommand {\semiF} {semi-Froeshl\'{e} } \newcommand{\Section}[1]{\setcounter{equation}{0}\section{#1}} \newtheorem{theorem}{Theorem} \newtheorem{corollary} {Corollary} \newcommand{\xvec} {{\bf x}} \newcommand {\tvec} {\mbox{\boldmath $\theta$}} \newcommand {\wvec} {\mbox{\boldmath $\omega$}} \newcommand {\chivec} {\mbox{\boldmath $\chi$}} \newcommand {\muvec} {\mbox{\boldmath $\mu$}} \newcommand {\xivec} {\mbox{\boldmath $\xi$}} \newcommand {\nvec} {{\bf n}} \newcommand {\mvec} {{\bf m}} \newcommand {\uvec} {{\bf u}} \newcommand {\bvec} {{\bf b}} \newcommand {\cvec} {{\bf c}} \newcommand {\gvec} {{\bf g}} \newcommand {\zvec} {{\bf z}} \newcommand {\Bvec} {{\bf B}} \newcommand {\upvec} {\mbox{\boldmath $\upsilon$}} \newcommand {\figdenominator} {(1)~} \newcommand {\figcontour} {(2)~} \newcommand {\figslopes} {(3)~} \newcommand {\figkplota} {(4)~} \newcommand {\figkplotb} {(5)~} \newcommand {\figkplotc} {(6)~} \newcommand {\figlogplot} {(7)~} \newcommand {\figbns} {(8)~} \newcommand {\figtda} {(9)~} \newcommand {\figtdb} {(10)~} \documentstyle[12pt]{article} \input{amssym.def} \input{amssym} \title{Breakup of Invariant Tori for the Four Dimensional Semi-Standard Map} \author{Erik M. Bollt \and James D. Meiss \\ \and Program in Applied Mathematics\\ Box 526, University of Colorado\\ Boulder CO 80309, USA} \date{October 1, 1992} \begin{document} \setlength{\baselineskip}{20pt} \maketitle \begin{abstract} We compute the domain of existence of two-dimensional invariant tori with fixed frequency vectors for a four-dimensional, complex, symplectic map. The map is a generalization of the semi-standard map studied by Greene and Percival; it has three parameters, $a_1$ and $a_2$ representing the strength of the kicks in each degree of freedom, and $\epsilon$, the coupling. The domain of existence of a torus in $(a_1,a_2)$ is shown to be complete and log-convex for fixed $k = \epsilon/a_1 a_2$. Explicit bounds on the domain for fixed $k$ are obtained. Numerical results show that quadratic irrationals can be more robust than the cubic irrational, ``the spiral mean." \end{abstract} \setlength{\baselineskip}{20pt} \newpage %\input{inc} %\input{psfig} %\include{intro} %\include{coupling} %\include{incommensurate} %\include{recursion} %\include{holomorphy} %\include{numerics} %\include{conclusions} %\include {fig} %\bibliography{JDM} %\end{document} \section{Introduction} Stability of motion in Hamiltonian systems and symplectic mappings is of great interest in many physical situations such as plasma and accelerator confinement and stellar and planetary dynamics; an understanding of stability is also of intrinsic theoretical interest. The primary stability result is the KAM theorem which asserts that most of the invariant tori of a nonlinear integrable Hamiltonian survive upon a {\it small}, smooth perturbation \cite{Arnold}. The robust tori, according to the theorem, are those which have sufficiently incommensurate frequency vectors (they satisfy a Diophantine condition, see \S{sec:incommensurate}). As a practical result, however, the KAM theorem has several drawbacks. The first is that estimates of the perturbation size for the destruction of tori are typically extremely small: much smaller than the size indicated by numerical computations for specific perturbations on specific tori (of course the theorem guarantees the survival of any Diophantine torus for any small enough perturbation). The second is that the theorem guarantees stability only for systems of two degrees of freedom since the invariant tori have half the dimension of phase space. None-the-less, computations indicate that while a system of three degrees of freedom may not be rigorously stable, it exhibits a ``practical stability" since orbits appear to remain trapped near invariant tori for extremely long periods. To some extent this is addressed by the Nekhoroshev theorem \cite{Nekhoroshev}, though this theorem requires extremely small perturbation sizes as well. For the case of two degrees of freedom, or equivalently area preserving mappings, much progress has been made in determining the existence of invariant tori. Three basic techniques have been used. The first is to examine the stability of a sequence of periodic orbits whose frequencies limit on the irrational frequency of interest---this gives rise to the {\it residue criterion} \cite{Greene,MacKay82}. It yields extremely accurate values for the parameters at which invariant circles are destroyed and can be made rigorous \cite{MacKay91}. The second method is a nonexistence criterion for twist maps, called converse KAM theory \cite{Mather,MacKay-Percival}. The final technique is numerical computation (in some cases using interval arithmetic) of the conjugacy to pure rotation \cite{Percival,DeLaLave,Berretti-Marmi,Marmi-Stark}. These methods can give accurate nonrigorous values for the critical parameter for essentially arbitrary Diophantine frequencies, and can also give reasonable rigorous values. Though many have attempted to generalize these techniques to Hamiltonian systems with more than two degrees of freedom, or equivalently, symplectic maps of four or more dimensions, there has been limited success in determining the existence of invariant tori. Periodic orbit approximations to invariant tori have been obtained \cite{Kook-Meiss,Muldoon}, and computations reveal that the stability domains of periodic orbits limiting on an incommensurate frequency vector may be converging for a volume preserving example \cite{Mao-Helleman1}. However the existence of the limit is difficult to prove \cite{Bernstein-Katok}, primarily because the ordering property of periodic orbits on the circle no longer applies on the torus. Converse KAM theory can be generalized to higher dimensions \cite{MacKay-Meiss-Stark}, though in this case one must assume that the tori are Lagrangian graphs. One of the fundamental problems in these studies is number theoretic: there is no satisfactory generalization of the continued fraction theory to simultaneous approximation of several irrationals (perhaps the most promising is that of Brentjes \cite{Brentjes}). In the case of the residue criterion, it is the best approximants (convergents of a continued fraction) whose properties converge to those of the invariant circle. Furthermore, quadratic irrationals play a large role in these studies because their continued fraction expansions are eventually periodic (these give rise to self-similar structures). Finally, the most robust tori appear to correspond to the class of quadratic irrationals known as the {\it noble} numbers; these have a continued fraction expansion with a tail of all one's. Roughly speaking, the explanation for this is that the noble numbers are the most difficult to approximate in the sense of Diophantine. The generalization of this class to higher dimensions is unknown. There has been some speculation that for four-dimensional mappings, cubic irrationals will replace the quadratics. One reason for this is that a periodic approximation scheme based on a Farey tree construction necessarily leads to a frequency which is the eigenvalue of a $3\times3$ matrix, and is therefore cubic \cite{Guckenheimer,Kim-Ostlund}. However, even in this case it has been difficult to determine if there is self-similar behavior near break-up \cite{Mao-Helleman1}, and there is no evidence that cubic irrationals are more robust than others. In this paper we study the four-dimensional, complex, symplectic map corresponding to the coupling of two semi-standard maps, as introduced in \cite{Greene-Percival}. This map is the complex version of a mapping introduced by Froeshl\'{e} \cite{Froeshle,Kook-Meiss}---we call it the \semiF map. We generalize the method of Percival and Greene \cite{Greene-Percival} to this case and find recursion formulae for the Fourier coefficients of an invariant two torus with a fixed frequency vector in \S{sec:recursion}. Existence of such a torus for small enough parameter values is guaranteed providing the frequency vector satisfies a Diophantine condition; we discuss this in \S{sec:incommensurate}. Because of the simple structure of the Fourier series for the \semiF map, we are able to apply some results from the theory of holomorphic functions of several complex variables in \S{sec:holomorphy}, and show that the domain of convergence of the Fourier series has a particular form; it is complete and log-convex. Finally in \S{sec:numerics} we compute these convergence domains for several example frequency vectors, including quadratic and cubic irrationals. \section{Incommensurate Frequencies} \label{sec:incommensurate} The convergence of the Fourier series for the semi-standard map has been studied extensively in \cite{Greene-Percival,Percival,Marmi-Stark}. In particular, rather sophisticated techniques for computing convergence of this series were developed in \cite{Percival}; these give accurate results for quite general frequencies. In general one determines a parameter interval $ |a| < a^{ss}(\omega)$ for which there is an analytic invariant circle with frequency $\omega$. Here $a^{ss}$, the critical function, is zero for every rational value and exhibits a maximum for \begin{equation} \label{golden} \omega = \gamma \equiv \frac{1+\sqrt{5}}{2} \;. \end{equation} The critical function appears to have a local maximum at each of the {\it noble frequencies}: those equivalent to $\gamma$ under a modular transformation, or equivalently which have a continued fraction expansion whose elements are all $1$ beyond some level. These results also apply to the \semiF map when $\epsilon =0$. Thus an invariant torus of frequency $\wvec = (\omega_1,\omega_2)$ exists within the rectangle $\{(a_1,a_2,\epsilon) : |a_1| < a^{ss}(\omega_1), |a_2| < a^{ss}(\omega_2), \epsilon = 0 \}$. Furthermore, since the \semiF map is an analytic perturbation of a twist map, KAM theory implies that for sufficiently small values of the three parameters $(a_1,a_2,\epsilon)$ there exists an invariant torus analytically conjugate to the rotation $\tvec \mapsto \tvec + 2\pi \wvec$ providing the frequency vector satisfies a Diophantine condition. For $d$-dimensions, the set of Diophantine vectors ${\cal D}_\mu$ consists of those $\wvec \in \Real^d$ for which there exists a $C > 0$ such that for all $({\bf p},q) \in \Integer^{d+1}$ \begin{equation} \label{Diophantine} |{\bf p} \cdot \wvec - q| \ge \frac{C}{||{\bf p}||^{\mu}} \;, \end{equation} where $||{\bf p}|| = \max (|p_1|,...,|p_d|)$. It is easy to see that if $\mu > d$, the measure of ${\cal D}_\mu$ approaches one as $C \rightarrow 0$; however, the measure of ${\cal D}_d$ is zero. Certainly if $\wvec \in {\cal D}_\mu$ then it is {\it incommensurate}, that is $1$ and $\omega_1, ...,\omega_d$ are linearly independent over the rationals. For $d>1$ one must distinguish between commensurate vectors and {\it resonant} vectors. While the former satisfy {\it some} rational relation ${\bf p} \cdot \wvec = q$, the latter have {\it all} components rational and correspond to periodic orbits. A straightforward generalization of Greene's method \cite{Greene} to higher dimensions would use resonant vectors, e.g. \cite{Guckenheimer}, instead of commensurate vectors. However, in KAM theory it is commensurabilities which cause the problems, not just resonances. Though there exist many Diophantine vectors, a result of Minkowski implies that every $\wvec$ can be closely approximated in a certain sense \cite{Cassels}: \begin{theorem} \label{minkowski} For any $\wvec \in \Real^d$ there are infinitely many integer vectors $({\bf p},q)$ such that when $K=1$ \begin {equation} |\bf p\rm \cdot \wvec -q|< {{K} \over {||\bf p\rm ||}^{d}} \;. \end{equation} If $d=1$ then $K$ can be replaced by $ 1/ \surd 5$ but nothing smaller. \end{theorem} To our knowledge, the minimal value of $K$ for $d>1$ is not known. One class of frequency vectors which are Diophantine are those constructed from algebraic irrationals \cite{Cassels}: \begin{theorem} \label{algebraic} If the components of $\wvec$ are incommensurate and elements of a real algebraic field of degree $d+1$, then $\wvec \in {\cal D}_d$. \end{theorem} Recall that an algebraic field generated by $\xi \in \Real$ of degree $n$ is defined as the set of numbers of the form $$R(\xi) = {P(\xi) \over Q(\xi)}$$ where P and Q are polynomials of degree $n$ with integer coefficients. One would expect that a frequency vector $\wvec$ which is more incommensurate, in the sense of having a larger Diophantine constant $C$ and smaller exponent $\mu$ would tend to persist for higher perturbations. This is numerically verified for the standard and semi-standard maps where the noble numbers give local maxima of $a^{ss}$, and are also the ``most" irrational in the sense of Diophantine. Unfortunately, to our knowledge, there are no results in the theory of simultaneous approximations which determine a class of frequency vectors analogous to the noble numbers. Indeed one of the main reasons for our numerical investigation is to attempt to develop a technique for determining this class. We will choose several simple frequency vectors as examples for our study. In addition to the golden mean, we will use the quadratic irrationals \begin{equation} \label{silver} \sigma \equiv \sqrt{2} = [2,2,2,2....] \equiv [2^{\infty}] \;; \ \ \zeta \equiv \frac{1+\sqrt{2}}{5+4\sqrt{2}} = [ 0,4,2^{\infty}]\;. \end{equation} The expressions on the right above give the continued fraction expansions. Setting $\wvec = (\gamma,\sigma)$ or $(\gamma, \zeta)$ yields two incommensurate frequency vectors since $\sqrt{{5 \over 2}}$ is irrational. Furthermore by Theorem \ref{algebraic}, both of these vectors are in ${\cal D}_2$, since they are elements of the algebraic field of degree three generated by $\xi = \sqrt{2} + \sqrt{5}$. This is easy to see, since any cubic polynomial in $\xi$ has the form $P(\xi)= a + b\sqrt{2} + c\sqrt{5} + d\sqrt{10}$ for $a,b,c,d \in \Integer$. Thus $\gamma$, $\sigma$, and $\zeta$ are all in $R(\xi)$. Finally we consider a cubic irrational, the real solution of \begin{equation} \label{eq:spiral} \begin{array}{lcl} \tau^3 &=& \tau+1 \\ \tau &\simeq& 1.32471795724474602596090885447809734 \\ &\simeq& [1,3,12,1,1,3,2,3,2,4,2,141,80,2,5,1,2,8,2,1,1,3,1,...] \;. \end{array} \end{equation} This so called ``spiral mean'' frequency, was introduced in \cite{Kim-Ostlund} as a possible analogue of the golden mean since in the Ostlund-Kim version of the Farey tree, $\tau$ has a simple periodic construction. The number $\tau$ is Diophantine since according to a theorem of Roth, every algebraic irrational is in ${\cal D}_{1+\delta} \; \forall \delta > 0$. Thus the critical function $a^{ss}(\tau) \ne 0$; however, determining its value is difficult because the continued fraction elements appear unbounded \cite{Mao-Helleman2}. None-the-less, for the four-dimensional case we will study the vector $( \tau, \tau^2)$ which is in the cubic field generated by $\tau$, and so an element of ${\cal D}_2$. Furthermore $\tau$ is the smallest of the ``{\it PV} numbers," which implies that the rational vectors on the Farey sequence approaching $(\tau, \tau^2)$ converge more slowly than any other algebraic pair \cite{Kim-Ostlund}. As we will see in \S{sec:numerics}, the frequencies enter the Fourier expansion for $\xvec (\tvec )$ solely in terms of the small denominators \begin{equation} \label{denom} D_\nvec=4\sin^2( \pi \nvec \cdot \wvec) \;; \end{equation} that is, the ${\bf n}^{th}$ Fourier coefficient of $\xvec(\tvec)$ is divided by $D_\nvec$. For Diophantine frequency vectors, $D_\nvec$ is bounded from below; in fact \eq{Diophantine} implies that if $\wvec \in{ \cal D}_{2}$, then $1/D_\nvec < {\cal O}(||\nvec||^4)$. Unfortunately, there is no theory analogous to the continued fraction theory which provides the values of $\nvec$ for which there are large peaks in $1/D_\nvec$. In fig. \figdenominator we show a plot of the values of $\nvec$ for which $(D_\nvec(\gamma,\sigma) ||\nvec||^4)^{-1} > 1.0\times10^{-2}$ and $5.0\times10^{-5}$. As this figure shows, these peaks are quite isolated and rare. Thus, following the results for the semi-standard map, one would expect the Fourier coefficients to have similar isolated peaks, and for the convergence determination of the Fourier series to be quite delicate. However, as we will see in \S{sec:numerics}, this is fortuitously not the case. \section{Coupling of Two Semi-standard Maps} \label{sec:coupling} The {\it semi-standard}, area preserving map was introduced by Greene and Percival \cite{Greene-Percival} as a numerically simpler model than the standard map for the investigation of the analytic properties of invariant circles. In Lagrangian form, the semi-standard map takes $\{x_{t-1},x_{t}\} \mapsto \{x_t,x_{t+1}\}$ and is defined by \begin{equation} \label {semistd} \delta^2 x_t \equiv x_{t+1}-2x_t+x_{t-1}=i a e^{ix_t} \;; \end{equation} this is a map on $\Complex^2$. The notation $\delta^2$ is reminiscent of the second derivative operator. In this paper we study a four dimensional generalization, analogous to the map introduced by Froeshl\'{e} \cite{Froeshle,Kook-Meiss}. Letting $\xvec_t \in \Complex ^2$, the \semiF map is \begin{equation} \delta^2\xvec_t \equiv \xvec_{t+1}-2\xvec_t+\xvec_{t-1}= {\bf F}(\xvec_t) \;, \label{semiF} \end{equation} where \begin{equation}\label{force} {\bf F}(\xvec) \equiv i \left( \begin{array}{ccc} a_1e^{ix^{(1)}}&+& \epsilon e^{ix^{(1)}+ix^{(2)}} \\ a_2e^{ix^{(2)}}&+& \epsilon e^{ix^{(1)}+ix^{(2)}} \end{array} \right) \;. \end{equation} There are three parameters, the strength of the kicks for each component semi-standard map ($a_1, a_2$) and $\epsilon$, the strength of the coupling of the two maps. \Eq{semiF} is symplectic since ${\bf F}$ is the gradient of a scalar potential (see for example \cite{Kook-Meiss}). We are looking for solutions $\xvec_t$ of \eq{semiF} which lie on an invariant two-torus homotopic to the trivial torus defined by the momentum ${\bf y}_t \equiv \xvec_t -\xvec_{t-1}$ being constant. In fact we demand that this torus be analytically conjugate to a uniform rotation on the angle variable $\tvec$ with a given frequency vector $\wvec$. These tori include those found by KAM theory. The conjugacy is represented by the following commuting diagram \begin{equation}\label{eq:3} \begin{array}{clcr} \xvec_t & \longrightarrow & \xvec_{t+1} \\ \downarrow & \ & \downarrow \\ \xvec(\tvec) & \longrightarrow & \xvec(\tvec+2 \pi\wvec) \; . \end{array} \end{equation} Thus, for a given $\wvec$, an invariant torus for \eq{semiF} is given by \begin{equation} \label{eq:4} \xvec_t = \xvec(\tvec+2\pi \wvec t) \;, \end{equation} for $\tvec \in \Torus^2$. The homotopy condition implies that \begin{equation} \label{eq:5} \xvec(\tvec+2\pi\mvec)=\xvec(\tvec)+2\pi\mvec \ \ \forall \mvec \in \Integer ^2 \;, \end{equation} thus $\xvec(\tvec)$ is coperiodic with $\tvec$; \begin{equation} \label{coperiodic} \xvec(\tvec)=\tvec+\chivec(\tvec) \end{equation} where $\chivec(\tvec)$ is doubly $2\pi$ periodic. If we suppose that $\xvec$ is analytic, it can be expanded in a Fourier series \begin{equation} \label{Fourier} \xvec(\tvec)=\tvec+\sum_{\nvec \in \Integer^{2}} \chivec_\nvec e^{i \nvec \cdot \tvec} \end{equation} Inserting \eq{eq:4} into \eq{semiF} yields the Percival form of the mapping \begin{equation}\label{Froeshle} \delta^2\xvec(\tvec) \equiv \xvec(\tvec+2 \pi \wvec)- 2\xvec(\tvec)+\xvec(\tvec- 2 \pi \wvec)= {\bf F}(\xvec(\tvec)) \; . \end{equation} Inserting the series (\ref{Fourier}) into \eq{Froeshle} will yield equations determining the Fourier coefficients $\chivec_\nvec$; these will be obtained in \S{sec:recursion}. \section{Recursion Relation} \label{sec:recursion} In this section we will derive recursion relations for the Fourier coefficients of ${\bf x}(\tvec)$, the solution to \eq{Froeshle}. For the semi-standard map it was possible to find a solution analytic in the upper half $\theta$ plane. In this case only the positive Fourier coefficients are nonzero. This is one advantage over the series for real mappings where all the Fourier coefficients must be considered. In the case at hand, since the force, \eq{force}, has only positive imaginary exponentials, we can also find solutions analytic in the domain $\{(\theta_1,\theta_2) : \rm{Im}{(\theta_1)} \geq 0 , \rm{Im}{(\theta_2)}\geq 0 \}$ , so that only positive Fourier coefficients are needed It is convenient to define $\uvec \in \Complex^2$ as \begin{equation} \label{complex} \uvec= \left( \begin{array}{c} a_1 e^{i\theta_1} \\ a_2 e^{i \theta_2} \end{array} \right) = \left( \begin{array}{cc} u_1 \\ u_2 \end{array} \right) \;. \end{equation} The advantage of this definition, is that the parameters $a_1$ and $a_2$ will not appear in any of our expansions. Further, using \eq{coperiodic} we define \begin{equation} \label{eq:12} {\bf g}(\uvec)=i({\bf x}(\tvec)-\tvec)=i \chivec (\tvec) \;. \end{equation} Since by ansatz, only the positive coefficients will be needed in the Fourier expansion of $\chivec (\tvec)$, ${\bf g}(\uvec)$ has a Taylor expansion \begin{equation}\label{taylor} {\bf g}(\uvec)=\sum_{\nvec\in \Natural^2} \bvec_\nvec \uvec^\nvec \equiv \sum_{n_1=0}^{\infty} \sum_{n_2=0}^{\infty} \left( \begin{array}{c} b_{(n_1,n_2)}^{(1)} \\ b_{(n_1,n_2)}^{(2)} \end{array} \right) u_1^{n_1} u_2^{n_2} \;. \end{equation} Here we use standard multi-index notation for the vector exponentiation: while $\uvec \in \Complex^2$ and $\nvec \in \Natural^2$, $\uvec^\nvec \equiv u_1^{n_1} u_2^{n_2} \in \Complex$. In addition to the expansion of $\bf g$, we will need the expansions of its exponential as well: \begin{equation} \label{exponential} e^{g_i(\uvec)}=\sum_{\nvec \in \Natural^2} c_\nvec^{(i)} \uvec^\nvec \end{equation} where $i = 1,2$. In terms of the new variables, the map \eq{Froeshle} takes the form \begin{equation} \label{eq:16} \delta^2 {\bf g}(\uvec)= - \left( \begin{array}{c} u_1 e^{g_1(\uvec)} \\ u_2 e^{g_2(\uvec)} \end{array} \right) - k\left( \begin{array}{c} u_1u_2 e^{g_1(\uvec)+g_2(\uvec)}\\ u_1u_2 e^{g_1(\uvec)+g_2(\uvec)} \end{array} \right) \end{equation} where \begin{equation}\label{k-define} k=\frac{\epsilon}{a_1a_2} \end{equation} is the coupling parameter. Note that these equations depend upon the three parameters $a_1,a_2$ and $\epsilon$ solely through $k$. Substituting \eqs{taylor}-(\ref{exponential}) into \eq{eq:16} and noting that for a term in the Fourier series the operator $\delta^2$ becomes $-D_\nvec$, as defined by \eq{denom}, yields the recursion relation for $\bvec_\nvec$ : \begin{equation} \label{b-recursion} D_\nvec \bvec_\nvec= \left( \begin{array}{c} c_{\nvec-(1,0)}^{(1)} \\ c_{\nvec-(0,1)}^{(2)} \end{array} \right) + k \sum_{\mvec=(0,0)}^{\nvec-(1,1)} \left( \begin{array}{c} c_\mvec^{(1)}c_{\nvec-\mvec-(1,1)}^{(2)}\\ c_\mvec^{(1)}c_{\nvec-\mvec-(1,1)}^{(2)} \end{array} \right) \;. \end{equation} If $\wvec$ is incommensurate, then $D_\nvec$ is nonzero, so that \eq{b-recursion} defines $\bvec_\nvec$. In fact $\bvec_\nvec$ is a convolution sum of $\{\cvec_\mvec\}$ for those $\mvec \prec \nvec$ where we define the partial order $\prec$ on integer vectors by $\mvec \prec \nvec$ if $m_i \le n_i$, and $\mvec \ne \nvec$. A simple derivative identity allows us to find the $\cvec_\mvec$ coefficients. \begin{equation} \frac{d}{du_j}e^{g_i(\uvec)}=[\frac{d}{du_j}g_i(\uvec)]e^{g_i(\uvec)}, \end{equation} which upon substitution of \eqs{taylor}-(\ref{exponential}) yields \begin{equation} \label{c-recursion} n_jc_\nvec^{(i)}=\sum_{\mvec \ne{(0,0)}}^\nvec m_j b_\mvec^{(i)} c_{\nvec-\mvec}^{(i)} \;. \end{equation} Note that \eq{c-recursion} allows the two forms, j=1 or 2, for $\nvec$ off the axis (these are equivalent), but for $\nvec$ on the axis, only one is valid because of a required division by a zero value of $n_j$. Examining \eq{c-recursion} reveals that $\cvec_\nvec$ is a function of strictly previous $\cvec_\mvec$, but up to current $\bvec_\nvec$, therefore the process must be started by generating $\bvec_\nvec$. Since the choice of initial phase $\tvec$ is arbitrary, we can set \begin{equation}\label{b0} \bvec_{\bf 0} = {\bf 0} \;. \end{equation} Examination of the mapping \eq{eq:16} yields in addition \begin{equation} \label{boundaries} b^{(2)}_{(n_1,0)}= 0 \;; \ \ b^{(1)}_{(0,n_2)} = 0 \;. \end{equation} Similarly, \eq{b0} and \eq{exponential} imply that $\cvec_{\bf 0} = {\bf 1}$, and \eq{c-recursion} yields \begin{equation} c^{(2)}_{(n_1,0)}= 0 \; ; \ \ c^{(1)}_{(0,n_2)}= 0 \; . \end{equation} Finally, the recursion (\ref{b-recursion}) implies that the values $b^{(1)}_{(n_1,0)}$ and $b^{(2)}_{(0,n_2)}$ are identical to those for the semi-standard map with frequencies $\omega_1$ and $\omega_2$, respectively. This completes the recursion algorithm which allows $\bvec_\nvec$ to be built as an explicit function of previous $\bvec_\nvec$ and $\cvec_\nvec$ coefficients. Note that if $k>0$ then $\bvec_\nvec$ is positive and real, a big advantage in their computation. Since equation \eq{taylor} actually represents two series, one in each component of the vector ${\bf g}$, the domain of convergence of ${\bf g}(\uvec)$ is the intersection of the domains of convergence of each component's series. \section{The Domain of Convergence} \label{sec:holomorphy} In this section, we review some relevant results on the domain of convergence for power series in several complex variables. Let $\zvec= (z_1,...z_d) \in \Complex^d$, and for $\mvec \in \Natural^d$, define ${\bf z}^m = z_1^{m_1}z_2^{m_2}....z_d^{m_d} \in \Complex$. We consider a power series, \begin{equation} S=\sum_{\mvec \in \Natural^d}^{} {b}_{\mvec}\zvec^\mvec \;, \end{equation} similar to the series obtained in the previous section. We denote the radii by $r_j = |z_j|$. The projection onto the radius space is denoted $\Pi: \Pi(\zvec) = (r_1, r_2,....r_d)$. The subset $\Complex^{d*} = \{\zvec : z_j \ne 0\}$ excludes points for which any component of $\zvec$ is zero. Several types of subsets of $\Complex^d$ are of interest. The {\it domain of convergence} of a series is the interior of the set of points for which it converges absolutely. A {\it polydisk} is the appropriate generalization of a disk: $P(a) = \{\zvec : |z_j| < |a_j|,\; j = 1,...,d\}$. A {\it Reinhardt domain} is a domain $R$ such that $R = \Pi^{-1}(\Pi(R))$; that is, if it contains a point with radii $r_j$, then it must contain every point with those same radii, regardless of phases. Reinhardt domains are conveniently pictured in the radius space $\Pi(\Complex^d) = \Real^d$. A Reinhardt domain is {\it complete} if for every $\zvec \in R$, the polydisk $P(\zvec) \subset R$; thus a complete domain contains all points with smaller radii. Finally a domain $D$ is {\it log-convex} if the set \begin{equation} \log(\Pi(D)) \equiv \{(\log(r_1), \log(r_2),....\log(r_d)): z \in \Complex^{d*} \cap D\} \end{equation} is a convex subset of $\Real^d$. We will use the following theorem \cite{Range,Kaup}: \begin{theorem} \label{th:convergence} If $S$ converges for all orderings of its terms at a point $\zvec$ then it converges absolutely to a holomorphic function. The domain of convergence, $D$, of $S$ is the interior of the set for which $|b_\mvec \zvec^\mvec|$ is bounded. Furthermore $D$ is a log-convex, complete Reinhardt domain. Conversely, if $|b_\mvec \zvec^\mvec|$ is unbounded then there is an ordering of the terms in $S$ for which it diverges. \end{theorem} The proof of this theorem is straightforward. Its most unusual aspect is that the domain of convergence is log-convex, which we will discuss in more detail. Suppose $\zvec, \xvec \in \Complex^{d*} \cap D$. Then for $\alpha + \beta = 1$ let $\uvec$ be any point in $\Complex^{d*}$ such that \begin{equation} \Pi(\uvec) = (r_1^\alpha s_1^\beta, ...,r_d^\alpha s_d^\beta) \;, \end{equation} where $r_j$ and $s_j$ are the radii of $\zvec$ and $\xvec$, respectively. Then, since $S$ converges at both $\zvec$ and $\xvec$, $B = \sup (|b_\mvec \zvec^\mvec|, |b_\mvec \xvec^\mvec|)$ exists, and \begin{equation} |b_\mvec \uvec ^\mvec| = |b_\mvec| \prod_{i=1}^{d} r_i^{m_i \alpha} s_i^{m_i \beta} \le B^{\alpha+\beta} = B \end{equation} is bounded as well. Thus $S$ converges at $\uvec$. Now since \begin{equation} \log(\Pi(\uvec)) = \alpha \log(\Pi(\zvec)) + \beta \log(\Pi(\xvec)) \;, \end{equation} we have shown that $D$ is log-convex. The application of Theorem \ref{th:convergence} to our system is straightforward since the series \eq{taylor} has the desired form; it yields the interesting result \begin{corollary}\label{convergence2} For fixed $k$ defined by \eq{k-define}, an analytic invariant torus with Diophantine frequency $\wvec$ of the \semiF map exists in a parameter domain in $(a_1,a_2)$ which is complete and log-convex. \end{corollary} In particular for fixed $k$, completeness implies that the domain of convergence is simply connected, and its boundary projected onto the radius space can be expressed as a graph of a function $r_1(r_2)$ or $r_2(r_1)$. As we will see in the next section, the calculation of these domains is possible with reasonable accuracy using the requirement that the terms in the series must be bounded. \section{Numerical Results} \label{sec:numerics} Determination of the sequence of $\{\bvec_\mvec\}$ of Fourier coefficients of $\chivec(\tvec)$ using the recursion algorithm of \S{sec:recursion} is straightforward, since they are real and positive for $k \geq 0$. The next issue is to numerically find the domain of absolute convergence, which from \S{sec:holomorphy}, is the set of $\uvec \in \Complex^2$ for which $|\bvec_\mvec \uvec^\mvec|$ is bounded. We begin by noting that the series (\ref{taylor}) \begin{equation} \label{g-series} \gvec(\uvec)=\sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \bvec_{m,n} u_1^m u_2^n \end{equation} converges absolutely in the polydisk $P(\uvec)$ if the reordered series \begin{equation} \label{eq:slopesum} \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \bvec_{m,n} r_{1}^{m} r_{2}^{n} =\sum_{n=0}^{\infty} r_1^n \Bvec_n (s) \;. \end{equation} converges. Here we define the slope $s \equiv r_2/r_1$, and the diagonal coefficient, \begin{equation} \label{diagonal} \Bvec_n(s) \equiv \sum_{m=0}^{n} \bvec_{n-m,m}s^m \end{equation} which, most importantly, is expressed as a finite sum. It follows that the radius of convergence for the $i^{th}$ component is \begin{equation} \label{log-radius} \log(r_1^{(i)}(s)) = - \lim_{n \rightarrow \infty} \frac{ \log{ B_n^{(i)}(s)}}{n} \end{equation} for each fixed $s$. Since the domain of convergence is complete according to Theorem \ref{th:convergence}, $r_1(s)$ is a single valued function. To avoid numerical overflow for large $s^m$, we use \eq{diagonal} when $s \leq 1$, and use a corresponding formula with the slope defined as $r_1/r_2$ otherwise. Furthermore, we can take advantage of definition (\ref{diagonal}) by computing the coefficients $\bvec_\mvec$ in a triangular domain $m+n \le N$. This saves about a factor of ten in computing time over using the square domain. Computer memory constraints lead us to chose $N=255$ as our matrix dimension. The efficacy of this method depends upon estimating \eq{log-radius}, the asymptotic growth rate of the Fourier coefficients. We first consider $\wvec=(\gamma,\sigma)$ where the components were defined by \eqs{golden} and (\ref{silver}). Figure \figcontour is a logarithmic contour plot of $b_\mvec^{(1)}$. It can be seen that $\bvec_\mvec$ grows rapidly as $\mvec$ grows, indeed the maximal values of $b_\mvec$ in the figure are ${\cal O}(10^{150})$. This is a result of the recursion algorithm, which shows that $\bvec_\mvec$ is a combination of all the previous $\cvec_\mvec$. When $k \ge 0$, the $\cvec_\mvec$ coefficients are positive and $\cvec_\mvec > \cvec_\nvec$ for $\mvec \succ \nvec$. Whenever there is a near commensurability, $D_\mvec$ is small, and $\bvec_\mvec$ takes a sudden jump. This can be seen in the contour plot as a serrating of the contour lines. The neighboring coefficients for greater $\mvec$ are influenced by this jump, but the recursion algorithm serves to spread and dissipate the extra height. In other words, the coupling serves to dampen the commensurabilities. This partly accounts for the stepping up nature of the contour plot. Fortunately, the limit (\ref{log-radius}) is not as difficult to evaluate for this problem as it could be in general. It turns out that $\log{B_n}$ behaves quite linearly as a function of n; this can be seen clearly in fig. \figslopes. The small spikes visible on a given ``line" of $(\log [B_n(s)],n)$ are due to near commensurabilities. By contrast, the Fourier coefficients for the semi-standard map, $b_n$, depend only on the small denominator and the single previous coefficient $c_{n-1}$, (which is in turn implicitly a function of the coefficients $b_1 ,..., b_{n-1}$). Resonances are extremely important, and primary, secondary, and even tertiary prominences can be observed, so that the Fourier coefficients have an extremely spiked profile. In the four-dimensional case the coupling between the frequencies $\wvec$ appears to play the dominant role. Resonances gain and lose prominence in a delicate balancing of the coupling between frequencies, which can be seen as shadows of vertical lines in \figslopes. To determine the radius in \eq{log-radius} we performed a least squares fit of variable data sets. The top and bottom ends of the fit were allowed to float by $n=10$ points each, and the fit with the lowest residual was automatically chosen. This eliminates the problem of a given fit falling just above or below a resonance spike. RMS errors in the slope fit are typically $\sigma_c =0.003 \pm 0.001$, which leads us to expect at least 2 decimal accuracy in the $r_i$ values. Using the three frequency pairs, $\wvec=(\gamma,\sigma)$, $(\gamma,\zeta)$, and $(\tau,\tau^2)$, we generate the respective $(r_1,r_2) = (a_1,a_2)$ curves for various coupling constants k, and for each of the $B_n^{(i)}$ components. These domains of convergence are displayed in figs. \figkplota-\figkplotc. $B_n^{(1)}$ is represented as solid curves, and $B_n^{(2)}$ as dashed curves. Figure \figlogplot displays the $(a_1,a_2)$ curves for $\wvec=(\gamma,\sigma)$ on a log- log scale. This example shows that D is log-convex in accord with Theorem \ref{th:convergence}. The sharp bends seen in some of the curves are due to the regular spacing of angles on a grid, which the log scale makes especially prevalent near the axes. Here we will discuss the behavior of the curves for the first component, $B_n^{(1)}$ (solid curves). When $a_2\rightarrow 0$, $r_1$ must approach $a^{ss}(\omega_1)$ since the map \eq{eq:16} becomes uncoupled in this limit, and $B_n^{(1)}$ becomes the coefficients of the semi-standard map with frequency $\omega_1$. We call the $r_1$ axis the ``dominant axis" for $B_n^{(1)}$; similarly, the $r_2$ axis will be the dominant axis for $B_n^{(2)}$. This behavior can be seen in figs. \figkplota-\figkplotc as all the various curves intersect the dominant axis at $a^{ss}(\omega_1)$. For reference, Table 1 gives the critical values of the semi-standard map for the various frequencies. Note that the curves in fig. \figkplota and \figkplotc actually overestimate the correct values on the axis; for example in fig. \figkplota, the intersection with the $r_2$ axis occurs near 0.985, while Table 1 implies that the correct value is 0.966. This overestimate is due to the fact that we compute the coefficients only out to the $255^{th}$ Fourier coefficient, and that near the axes the spikes in the $B_n$ curves become more prominent (see further discussion of this below). For the semi-standard mapping more sophisticated fitting techniques (e.g. \cite{Percival}) are required for an accurate evaluation of the critical function. For our mapping we believe that, away from the axes, the radius curves are actually more accurate than this indicates. In fig. \figkplotb, the intersection with the $r_1$ axis appears to be much lower than the value $r_1 = 0.979$ given in Table 1; however these curves actually rise rapidly to the correct (actually overestimated) value as $r_2 \rightarrow 0$. It is interesting that in this case even though the values on axis are quite different, the convergence boundary has adjusted itself to be nearly square for small $k$. Finally, this rapid rise---approaching $a^{ss}(\omega_1)$ at a sharp angle, does not violate log-convexity, as required by Theorem \ref{th:convergence}. \begin{table} \label{table} \[\begin{array}{cl}\omega &{\rm a}^{\rm ss}\\ \gamma &\rm 0.979661\\ \sigma &\rm 0.966165\\ \zeta &\rm 0.833726\\ \tau &\rm 0.657\\ {\tau }^{2}&0.660\\ \gamma \rm +\sigma &\rm 0.09\\ \gamma \rm +\zeta &\rm 0.66\\ \tau \rm +{\tau }^{2}&0.33\end{array}\] \caption{Critical values for the semi-standard map.} \end{table} The figures also show that the solid curves limit to $a^{ss}(\omega_2)$ on the $r_2$ axis, which we call the ``subdominant" axis for $B_n^{(1)}$. This phenomena requires some explanation. When $\epsilon \equiv 0$, the boundary of domain of convergence for $B_n^{(1)}$ is $r_1^{(1)} = a^{ss}(\omega_1)$, independent of $r_2$; the numerical results for nonzero $\epsilon$, however, imply that $r_2$ limits to $a^{ss}(\omega_2)$ on the $r_2$ axis. This also occurs for the domain of convergence of the second component $B_n^{(2)}$: $r_1 \rightarrow a^{ss}(\omega_1)$ as $r_2 \rightarrow 0$. To explain this phenomena, consider for example the small slope limit of $B_n^{(2)}(s)$. \Eq{diagonal} implies that \begin{equation} {B}_{n}^{(2)}(s) = {b}_{n,1}^{(2)}\ s+\ {\cal O}({s}^{2}) \end{equation} where the $s^0$ term vanishes according to \eq{boundaries}. Using the recursion relation (\ref{b-recursion}) implies \begin{equation} \label{small-s} B_n^{(2)}(s) \simeq sk{\frac{D_{(n,0)}}{D_{(n,1)}}} b _{(n,0)}^{(1)} \end{equation} Thus using \eq{log-radius} the radius of convergence is \begin{equation} \log(r_1^{(2)}) \simeq \log \left[a^{ss}(\omega_1)\right]\ -\ \lim\limits_{n\ \rightarrow \ \infty} {\frac{1}{n}}\log\left({\frac{D_{(n,0)}}{D_{(n,1)}}}\right) \end{equation} The last limit in fact is zero, since by the Diophantine condition \eq{Diophantine}, the ratio of the denominators is bounded by ${\cal O}(n^4)$. A similar result holds for the first component of the mapping along the $r_2$ axis, so we have shown that \begin{equation} \label{subdom} \lim\limits_{s\ \rightarrow \rm \ 0}^{}\ {r}_{1}^{(2)}(s)\ =\ {a}^{ss}(\omega _1),\ \ \lim\limits_{s\ \rightarrow \rm \ \infty }^{}\ {r}_{2}^{(1)}(s)\ =\ {a}^{ss}(\omega _2) \end{equation} Furthermore \eq{subdom}, together with completeness implies that the the domain of convergence is bounded by the rectangle \begin{equation} a_1 \leq a^{ss}(\omega_1),\ \ \ a_2 \leq a^{ss}(\omega_2) \end{equation} In fact the figures show that as $k \rightarrow 0$ the domain of convergence approaches this rectangle. Our interpretation of this is that for small but nonzero k, the singularity corresponding to $r_2=a^{ss}(\omega_2)$ is still present, though weakened (the ``residue" of this singularity, limits to zero as $k \rightarrow 0$, but it is still present for any nonzero k). This causes a difficulty with our numerical scheme for finding $r_1(r_2)$ when $k$ is small; we we discuss this further below. Figure \figbns displays the coefficients $B_n^{(1)}$ and $B_n^{(2)}$ for $s= 10^{- 25}$ and $\wvec=(\gamma,\sigma)$. In the limit of small slope $B_n^{(1)} \approx b_{n,0}^{(1)}$ which are the Fourier coefficients for the semi-standard map \cite{Greene-Percival}. Thus the upper plot is indistinguishable from that for the semi-standard map. The lower half of fig. \figbns shows $B_n^{(2)}$ for small slope. The profile exhibits spikes and valleys corresponding to a complicated coupling between $\sigma$ resonances and the still important $\gamma$ resonances, as shown in \eq{small-s}. Furthermore, \eq{small-s} implies that the profile approaches a limiting form as $s \rightarrow 0$, even though the magnitude of $B_n^{(2)}$ approaches zero. Likewise, $B_n^{(2)}$ near the $r_2$, axis yields the semi-standard coefficients for $\omega_2 = \sigma$, while $B_n^{(1)}$ goes to zero, while similarly converging to a fixed profile. As the domain of convergence plots show, the rectangular domain for small $k$ contains the domain for any finite $k$. This follows from the completeness of the domain of convergence, and the fact that the curves limit to the semi-standard values on the axes. This fact can be used as an upper-bound when discussing the question of which torus is ``last." We also computed $r_1(r_2)$ curves for negative values of k. By the same argument as above, the negative k curves intersect the axis at $a^{ss}(\omega_1)$ and $a^{ss}(\omega_2)$. Otherwise the curves are qualitatively similar to those shown in figs. \figkplota-\figkplotc, so we omit the plots. Since the domain of convergence depends only on $k$, these curves provide the boundary of existence in four of the octants in $(a_1,a_2, \epsilon)$ space, the other four being determined by the positive $k$ results. As $k$ increases all of the boundaries in figs. \figkplota-\figkplotc become hyperbolic in shape. This can be seen most clearly in fig. \figkplotb, for $\wvec=(\gamma,\zeta)$. The large $k$ limit corresponds to \begin{equation} \epsilon \gg (a_1 , a_2). \end{equation} Taking this to the extreme, we set $a_1=a_2=0$, then \eqs{semiF}-(\ref{force}) have the form \begin{equation} \label{eq:hyp} \delta^2 \xvec=i\epsilon \left( \begin{array}{c} e^{ix_1+ix_2}\\ e^{ix_1+ix_2} \end{array} \right). \end{equation} Defining the new variables \begin{equation} \xi_1=x_1+x_2 \end{equation} \begin{equation} \xi_2=x_1-x_1, \end{equation} and adding and subtracting the components of \eq{eq:hyp} yields a new map. \begin{equation} \label{eq:hypmap} \begin{array}{l} \delta^2 \xi_1= 2i \epsilon e^{i \xi_1}\\ \delta^2 \xi_2=0. \end{array} \end{equation} Thus, there exists an invariant torus for $(\xi_1,\xi_2)$ up to some critical value \begin{equation} \label{eq:maxeps} 2 \epsilon=a^{ss}(\omega_1+\omega_2). \end{equation} Now \eq{eq:hyp} is approximately valid for small $a_1$ and $a_2$, so we expect that as $k \rightarrow \infty$, using $\epsilon = k a_1 a_2$, the fixed $k$ boundary will limit to \begin{equation} \label{hyperbola} r_1r_2 = \frac{a^{ss}(\omega_1+\omega_2)}{2k} \end{equation} which defines a hyperbola. Thus we have three analytic bounds on the domain of existence of a torus: \begin{equation} \label{bounds} a_1 < a^{ss}(\omega_1), \ \ \ a_2 < a^{ss}(\omega_2), \ \ \ \epsilon < 0.5 a^{ss}(\omega_1+\omega_2) ~~, \end{equation} though the last equation is not rigorously derived. As a confirmation, Table 1 shows that $a^{ss}(\gamma+\sigma)$ is much smaller than $a^{ss}(\gamma+\zeta)$ and $a^{ss}(\tau+\tau^2)$. Thus, \eq{hyperbola} predicts that the curves for $\wvec=(\gamma,\sigma)$ should become hyperbolae more quickly than for other $\wvec$ curves, as we do in fact observe. As mentioned earlier, the scheme (\ref{eq:slopesum})-(\ref{log-radius}) for finding $r_1(s)$ has numerical problems when $k \ll 1$. For such small k, the singularity on one axis is dominant over the singularity on the other axis. To illustrate the problem, consider a simple example which has a similar imbalance in the prominence of its singularities. Let \begin{equation} \label{exampleseries} S(r_1,r_2)=\frac{\alpha}{\alpha-r_1} + \frac{\delta \beta}{\beta-r_2} = \sum_{m,n} b_{m,n}r_1^mr_2^n~~~. \end{equation} Here small values of $\delta$ simulate small values of $k$; however, for any nonzero $\delta$, the domain of convergence of this series is the rectangle $\{(r_1,r_2): r_1 < \alpha, r_2 < \beta\}$. We examine the behavior of equations (\ref{eq:slopesum})-(\ref{log-radius}) when applied to \eq{exampleseries} by a perturbation analysis near $s=0$. For a finite $n$, the algorithm gives an error in $r_1$ of \begin{equation} \label{error} \Delta r_1 \ \sim \ -\ {\frac{\delta \alpha }{\rm n}}\ {\left({{\frac{\alpha s}{\beta }}}\right)}^{n}~~~. \end{equation} Thus the method works well provided $s < \beta / \alpha$, but fails drastically for larger $s$. In our computations, the slope is never larger than one; we switch to the inverse of the slope when $s=1$. Thus, supposing $\beta < \alpha$ the method fails in a cone $\beta / \alpha < s < 1$. So for the Froeshl\'e mapping, we also expect that slopes within a similar cone will give bad results if k is too small. That this is true can be seen as a slight loss of convexity for the smallest values of k along the subdominant axis in figs. \figkplota- \figkplotc. In practice we are unable to lower $k$ below $10^{-5}$ in the computations. Finally, our $r_1(r_2)$ data can be displayed in terms of the coupling parameter $\epsilon$, instead of k. Figures \figkplotb and \figkplotc are converted via \eq{k-define} to the three dimensional graphs seen in figs. \figtda and \figtdb. Here we see in a new way the importance of the sum frequency $(\omega_1+\omega_2)$ through \eq{bounds}. Numerical overflow for large k prevents us from calculating the curves for $\epsilon$ too close to its maximum value. In many ways, it is these three dimensional plots which are most useful when deciding a partial order to determine the ``last invariant torus." One concept of ordering of the domains of convergence is to choose a directed curve in $(a_1,a_2, \epsilon)$ beginning at the origin. One could linearly order the domains of convergence in terms of the order of intersection of the domain boundaries with this curve. This motivates the following local definition of order. \bf{Curve based Order}: \it{An $\wvec$ torus persists longer than a $\muvec$ torus along a curve $\xivec(t)$ for which $\xivec(0) = \bf{0}$, if $\xivec(t)$ intersects the boundary of the domain of convergence of the $\muvec$ torus first.}\rm The simplest example of a parameterized family is a line emanating from the origin. Another example is a parabolic ray $a_1 = t, a_2 = st, \epsilon = k s t^2$ for fixed $s$ and $k$. Figures \figkplota-\figkplotc order domains in this sense. In general, one wants to do more than compare two surfaces using a single point from each surface, which is all a curve based order allows. In some sense, one may want to incorporate the information of the entire surface in a comparison. This motivates the definitions of the following global comparisons. \bf{Metric based Order}: \it{For a given metric, an $\wvec$-torus persists longer than a $\muvec$-torus if the boundary of the domain for $\wvec$ has a point farther from the origin than that for $\muvec$.}\rm This definition for ordering is limited in that it requires the choice of a metric. If one surface is completely contained inside another, then that torus is more persistent than the other according to any definition, since containment is a topological notion. Thus we define the \it{partial ordering} \bf{Topological Order}: \it{An $\wvec$ torus persists longer than a $\muvec$ torus if the domain for $\wvec$ contains that of $\muvec$.}\rm Of course, the surfaces for two different frequencies will intersect in general, and then the topological ordering does not apply. In our examples, the surface for $(\tau,\tau^2)$ is completely contained inside that of $(\gamma,\zeta)$, and therefore the $(\gamma,\zeta)$ torus is more persistent. The complete containment of the $(\gamma,\zeta)$ surface is partly due to the fact that each of $(a^{ss}(\gamma),a^{ss}(\zeta),a^{ss}(\gamma+\zeta))$ are greater than their counterparts $(a^{ss}(\tau),a^{ss}(\tau^2),a^{ss}(\tau+\tau^2))$. On the other hand in order to compare the $(\gamma,\sigma)$ and $(\gamma,\zeta)$ tori, note that though $a^{ss}(\sigma)) > a^{ss}(\zeta))$, $a^{ss}(\gamma+\sigma)< a^{ss}(\gamma+\zeta)$. Thus the surfaces must intersect, and therefore there can only be parameterized comparisons. \section{Conclusions} \label{sec:conclusions} We have determined the domain of existence of invariant two-tori analytically conjugate to a rotation for the \semiF mapping by expanding the conjugacy function in a Fourier series in the angle variables. The \semiF mapping has the advantage that two of the parameters can be eliminated in the Fourier series, so that the boundary of existence of the tori in all three parameters can be obtained with a single parameter sweep. We have studied the boundary of the domain for several frequency vectors, all of which are elements of a cubic algebraic field, and therefore satisfy Diophantine conditions. The boundary of these domains appears to be smooth; rather surprisingly, it appears smooth even when the parameters have opposite signs (i.e. negative values of $k$). We have shown that when projected on the parameters $(a_1,a_2)$ for fixed $k = \epsilon/a_1 a_2$, the boundary is log-convex and complete, and that as $k \rightarrow 0$ the domain limits to the rectangle corresponding to the domain for the uncoupled mappings. Furthermore, numerical results imply that the domain is bounded by the critical function for the sum frequency, as shown by \eq{eq:maxeps}. The methods and theorems of this paper are not restricted to the four-dimensional version of \eq{semiF}. They also apply to the 2d dimensional complex \semiF map, providing only that each occurrence of $x^{(j)}$ in an exponential, $\exp{(im x^{(j)})}$, in the force has the same sign. The main bottleneck is computing the Fourier coefficients recursively which involves a (d-1) degree iterated convolution sum, where d is the dimension. Computing the $m^d$ coefficients would take ${\cal O}(m^{2d})$ steps, making computer time a major problem in practice. In the same vein, more complicated forcing terms in \eq{force} could also be considered, but similar time constraints may be a problem. There are a number of open questions left by our study. 1) When the Fourier series does not converge, does there exist an invariant Cantor set for the mapping (a cantorus)? Results for twist mappings near the anti-integrable limit show the existence of cantori for all frequencies \cite{MacKay-Meiss}. What is the nature of the invariant set when the Fourier coefficients for $x^{(2)}$ converge, but those for $x^{(1)}$ do not, as seen especially in Fig. \figkplotc? One is tempted to think it is a Cantor set of circles. 3) Are all invariant tori for the \semiF mapping analytically conjugate to a rotation? Perhaps all tori with Diophantine frequency vectors? 4) Is there an extension of the converse KAM theory of \cite{MacKay-Meiss-Stark} to complex mappings? 2) Is there an extension of some of the results of Theorem \ref{th:convergence} to real valued four-dimensional mappings of some class? It is possible that such a map may also have a log-convex domain in the proper coordinates. 5) Can one use similar techniques to study the existence of invariant circles for a four-dimensional mapping? In \cite{MacKay-Meiss-Stark} it was suggested that circles may last longer than any tori. 6) Which class of frequency vectors correspond to the most persistent invariant tori? In this paper we compare several likely candidates, but do not present evidence that there are not more persistent tori. In searching for a particularly persistent torus, a first step might be to maximize the values of $a^{ss}(\omega_1), a^{ss}(\omega_2)$, and $a^{ss}(\omega_1+\omega_2)$. Which class of frequency vectors does this? Of course since denominators containing all $\mvec \cdot \wvec$ occur, the most persistent class of frequencies may be that with maximal Diophantine constant $C$ in \eq{Diophantine}. Since incommensurate algebraic frequency vectors form a field, any elements of such a field will have the same $C$. Moreover, since a degree three algebraic field has the minimal exponent $\mu$ in \eq{Diophantine}, it seems reasonable that it is such a field which will be most persistent. Of course the definition of persistence will depend on the choice of a partial ordering, and even then it is not clear how dependent upon the specific model the results would be. \section*{Acknowledgements} Partial support for this research was obtained from the NSF through grant DMS-9001103 and by the ONR through grant N00014-91-4037. We would aslo like to thank John Greene for helpful discussions regarding the convergence of the semi-standard map critical function. For figures, send to Erik Bollt Program in Applied Mathematics University of Colorado, Boulder 80309-0526 (303)492-4668 -4066FAX \begin{thebibliography}{10} \bibitem{Arnold} V. I. Arnold,{\it Mathematical Methods of Classical Mechanics}, (Springer, New York), 462pp. (1978). \bibitem {Berretti-Marmi} A. Berretti and S. Marmi, ``Standard Map at Complex Rotation Numbers: Creation of Natural Boundaries", \newblock {\it Phys. Rev. Let.} {\bf 68}, 1443-1446(1992). \bibitem{Bernstein-Katok} D. Bernstein and A. Katok, ``Birkhoff Periodic Orbits for Small Perturbations of Completely Integrable Hamiltonian Systems with Convex Hamiltonians", \newblock {\it Invent. Math.} {\bf 88}, 225(1987). \bibitem{Brentjes} A. J. Brentjes, ``Multidimensional Continued Fraction Algorithms", Ph.D. Thesis, \newblock Mathematisch Centrum, Amsterdam, (1981). \bibitem{Cassels} J. W. S. Cassels, {\it An Introduction to Diophantine Approximation}, (Cambridge University Press, Cambridge), 168 pp. (1965). \bibitem {DeLaLave} R. de la Llave and D. Rana, ``Accurate Strategies for Small Divisor Problems", \newblock {\it Bull.AMS} {\bf 22}, 85-90(1990). \bibitem{Froeshle} C. Froeschle, ``Numerical Study of a Four Dimensional Mapping", \newblock {\it Astron. and Astrophys} {\bf 16}, 172-189(1972). \bibitem{Froeshle2} C. Froeschle and J.-P. Scheidecker, ``Numerical Study of a Four Dimensional Mapping II", \newblock {\it Astron. and Astrophys.} {\bf 22}, 431-436(1973). \bibitem{Greene} J. M. Greene, ``A Method for Computing the Stochastic Transition", \newblock {\it J. Math. Phys.} {\bf 20}, 1183-1201(1979). \bibitem{Greene-Percival} J.~M. Greene and I.~C. Percival, ``Hamiltonian Maps in the Complex Plane", \newblock {\it Physica} {\bf 3D}, 530-548(1981). \bibitem{Guckenheimer} J. Guckenheimer, B. Hu and J. Rudnick, ``Quasiperiodic Transition to Chaos with Three Incommensurate Frequencies", \newblock unpublished, (1983). \bibitem{Kaup} L. Kaup and B. Kaup, {\it Holomorphic Functions of Several Variables}, (W. de Gruyter,Berlin), 349 pp. (1983). \bibitem{Kim-Ostlund} S. Kim and S. Ostlund, ``Simultaneous Rational Approximations in the Study of Dynamical Systems", \newblock {\it Phys. Rev. A} {\bf 34}, 3426-3434(1986). \bibitem{Kook-Meiss} H. T. Kook and J. D. Meiss, ``Periodic Orbits for Reversible, Symplectic Mappings", \newblock {\it Physica} {\bf 35D}, 65-86(1989). \bibitem{MacKay82} R. S. MacKay, ``A Renormalisation Approach to Invariant Circles in Area-Preserving Maps", \newblock {\it Physica} {\bf 7D}, 283-300(1983). \bibitem{MacKay91} R. S. MacKay, ``On Greene's Residue Criterion", \newblock University of Warwick, (1991). \bibitem{MacKay-Meiss} R. S. MacKay and J. D. Meiss, ``Cantori for Symplectic Maps near the Anti-integrable Limit", \newblock {\it Nonlinearity} {\bf 5}, 149-160(1992). \bibitem{MacKay-Meiss-Stark} R. S. MacKay, J. D. Meiss and J. Stark, ``Converse KAM Theory for Symplectic Twist Maps", \newblock {\it Nonlinearity} {\bf 2}, 555-570(1989). \bibitem{MacKay-Percival} R. S. MacKay and I. C. Percival, ``Converse KAM: Theory and Practice", \newblock {\it Comm. Math. Phys.} {\bf 98}, 469-512(1985). \bibitem {Malavasi-Marmi} M. Malavasi and S. Marmi, ``Critical KAM Circles and the Brjuno Function", \newblock {\it J. Phys. A} {\bf 22}, L563-L569(1989). \bibitem{Mao-Helleman1} J. M. Mao and R. H. G. Helleman, ``Existence of KAM Tori for Four-Dimensional Volume Preserving Maps", \newblock {\it N. Cimento} {\bf 104B}, 177-183(1989). \bibitem{Mao-Helleman2} J. M. Mao and R. H. G. Helleman, ``Break-up of KAM Tori of Cubic Irrational Winding Number", \newblock Univ. of Houston, (1989). \bibitem {Marmi-Stark} S. Marmi and J. Stark, ``On the Standard Map Critical Function," \newblock {\it Nonlinearity}, (in press). \bibitem{Mather} J. N. Mather, ``Non-Existence of Invariant Circles", \newblock {\it Erg. Th. Dyn. Sys.} {\bf 2}, 301-309(1984). \bibitem{Muldoon} M. Muldoon, ``Ghosts of Order on the Frontier of Chaos" \newblock Ph.D. Thesis, California Institute of Technology (1989). \bibitem{Nekhoroshev} N. N. Nekhoroshev, ``An Exponential Estimate for the Stability Time of Hamiltonian Systems Close to Integrable Ones", \newblock {\it Russ. Math. Surveys} {\bf 32:6}, 1-65(1977). G. Benettin, L. Galgani and A. Giorgilli, ``A Proof of Nekhoroshev's Theorem for the Stability Times in Nearly Integrable Hamiltonian Systems", \newblock {\it Celest. Mech.} {\bf 37}, 1-25(1985). \bibitem {Percival} I. C. Percival, ``Chaotic Boundary of a Hamiltonian Map", \newblock {\it Physica} {\bf 6D}, 67-77(1982). \bibitem{Range} R. M. Range, {\it Holomorphic Functions and Integral Representations in Several Complex Variables}, (Springer-Verlag, New York), 386 pp. (1986). \end{thebibliography} \end{document}