Testing the accuracy of the overlap criterion

Here we investigate the accuracy of the overlap criterion when applied to a simple near-integrable model in both its 2D and 3D version. To this end, we consider respectively, two and three quartic oscillators as the unperturbed system, and couple the degrees of freedom by a cubic, non-integrable perturbation. For both systems we compute the unperturbed resonances up to order O(\epsilon^2), and model each resonance by means of the pendulum approximation in order to estimate the theoretical critical value of the perturbation parameter for a global transition to chaos. We perform several surface of sections for the bidimensional case to derive an empirical value to be compared to our theoretical estimation, being both in good agreement. Also for the 3D case a numerical estimate is attained that we observe matches the critical value resulting from theoretical means. This confirms once again that reckoning resonances up to O(\epsilon^2) suffices in order the overlap criterion to work out. Keywords: {Chaos -- Resonances -- Theoretical and Numerical Methods}


I. INTRODUCTION
Though the stability problem of Hamiltonian systems has been almost completely elucidated by a rigourous sequence of theorems that build up the so called KAM theory (see for instance [1] and [2], together with the original references therein: [3], [4] and [5]), the application of the results of the KAM theory to a specific system is far from being an easy task. In fact, it turns out to be much simpler to take advantage of the heuristic Overlap Criterion, which seems to provide similar estimations to those resulting from the KAM theory.
The overlap criterion due to Chirikov (see [1]) has been largely used in many different fields, its probably most popular application being the study of inestabilities in the Solar System as well as in other planetary models (see for example [7], [8], [9]). In any case, since the widespread model for a resonance is the pendulum approximation, the overlap criterion applies directly to the intersection of their associated unperturbed separatrices (or heteroclinic intersections).
In his pioneer work on the standard map [1] Chirikov shows that the application of the overlap criterion to primary resonances overestimates the actual value of the critical parameter, K c , and only when high order resonances are considered, does the overlap criterion succeed in providing a more accurate value for K c . In fact, the author shows that on including the third harmonics resonances, the overlap criterion leads to K c ≈ 1, rather close to the empirical value.
In the present effort we address a similar analysis to that performed by Chirikov, but using a 2D and a 3D nearintegrable Hamiltonian systems, namely, two and three uncouppled quartic oscillators perturbed by a cubic term. The 3D version of this model has been studied in [10], [12], where the authors numerically investigate the global dynamical properties of the model and estimate the critical value ǫ c beyond which the system is globally chaotic, i.e. for which less than the 10% of the energy surface corresponds to invariant tori. By means of the overlap criterion we derive such a critical value for the perturbative parameter ǫ on considering not only primary but also high order resonances, and the perturbation Fourier series truncated at O(23 −2 ) in their coefficients. We then compare, for each case, the theoretical critical value with that obtained by numerical means.
The paper is organized as follows. The dynamical system under study in its 2D version is described in Section II and its relevant resonances at O(ǫ) are obtained in Section III, their widths being determined in Section IV. The resonances at O(ǫ 2 ) are provided in Section V, where an estimate of the critical value of the perturbative parameter is provided. For the sake of comparison, an empirical estimate of such a value is obtained in Section VI by recourse of performing several surfaces of section for the system. The 3D model is addressed in Section VII, whose resonances at order O(ǫ) and O(ǫ 2 ) are given in Sections VIII and IX respectively. Section X is devoted to the theoretical estimate of the critical parameter, which is shown to be in good agreement with the one given in [12]. A final discussion is provided in Section XI.

II. THE 2D DYNAMICAL SYSTEM
Here we will be concerned first with a 2D perturbed quartic oscillator. In cartesian coordinates the system is described by the following Hamiltonian (see [10]): where ǫ is a perturbative parameter that controls the strength of the perturbation. On setting ǫ = 0 we recover the integrable quartic oscillator's Hamiltonian (see [11], [10] and references therein), whose solutions are given by where we have used the following definitions: where K(k) denotes the complete elliptic integral, and the coeficients in the Fourier expansions (2) satisfy: The third equation in (3) reveals the dependence of the frequencies on the unperturbed energies, enabling us to get the functional relationship between the latter and the unperturbed action variables, namely, h i = AI i 4/3 , with A = (3β/2 √ 2) 4/3 . With this relation in mind, and taking into account that the angle variables are θ i ≡ ω i (h i )t, i = 1, 2, the complete Hamiltonian, in terms of the action-angle variables of the unperturbed Hamiltonian, can be recast as: where H 0 (I) = A(I 1 4/3 + I α nmk cos 2(n + m − 1)θ 1 ± (2k − 1)θ 2 + cos 2(n − m)θ 1 ± (2k − 1)θ 2 with α nmk ≡ α n α m α k , andV (I) ≡ 2 5/2 3β 4 I 2/3 1 I 1/3 2 , the ± sign meaning that both terms are included in the series.

III. RESONANCES AT O(ǫ)
A glance at the perturbation series given in equation (5) reveals that the number of resonant terms at first order in the perturbative parameter is unbounded, which is a drawback to take into account the width of every resonance at such an order.
However, the strong dependence of the Fourier amplitudes on (n + m + k), through the quantities α nmk ≈ 1/23 (n+m+k−3) , gives us a good hint on how to gather the O(ǫ) resonances and where to cut the series, our attempt being to keep terms only up to O(1/23 2 ).
All the possible combinations of n, m and k verifying that n + m + k ≤ 5, yield 24 different vectors which are listed in Table I, together with the number of times they appear in a term with coefficient α nmk of a given order in 1/23. Thus, N 0 denotes the number of times the vector appears with coefficient α nmk = α 3 1 , N 1 the number of times it arises with coefficient α 2 1 α 2 (≈ α 3 1 /23), and N 2 corresponds to either the coefficient α 1 α 2 2 or α 2 1 α 3 (which are approximately α 3 1 /23 2 ). From now on, these vectors will be appointed as harmonics at order O(ǫ, 1/23 2 ).  On applying the resonance condition m · ω = 0, with m ∈ Z 2 /{0}, to the unperturbed system, the following relation between the energies in each degree of freedom is obtained: which implies that the resonance structure in energy and action space consists of straight lines (with positive slope) given by and respectively. Moreover, eq. (6) indicates that m 1 m 2 ≤ 0, whence, those vectors having both components with the same sign must be discarded. Let us notice however, that not all of the remaining harmonics at order O(ǫ, 1/23 2 ) are actually resonant, as it will be discussed in the forthcoming section.

IV. WIDTH OF THE RESONANCES AT O(ǫ)
On computing the widths of resonances, the pendulum's approximation (see for instance [1], [13]) provides a suitable description whenever each resonance is assumed to be isolated from the rest.
Notice should be taken that, before proceeding to estimate the width of a given resonance, all the coefficients α nmk associated to the same trigonometric function are to be added together into a single one. Indeed, for each given vector m, we should define the coefficient: where n, m, and k are natural numbers that combine to form the vector m in any of the four ways displayed in eq. (5). Let us remark that α m ≈ α 3 1 (N 0 + N 1 /23 + N 2 /23 2 ). Inasmuch the pendulum's approximation has been applied to this single resonant term, the new (resonant) Hamiltonian turns out to be: with where the sum over repeated indexes should be understood. Let p r be the maximum variation of p 1 within the oscillation regime, then As a consequence of the simple pendulum dynamics, the maximum displacement of the unperturbed action variables depends on m and p r in the fashion: (∆I) r ≡ (I − I r ) max = p r m. Furthermore, the maximum displacement of the unperturbed energy is given by and it is the maximum amplitude attained in the oscillation of any of the unperturbed energies that measures the width of the resonance.
Let us now recall that in any 2D problem the energy conservation condition h = h 1 +h 2 , together with the resonance condition m 1 h 1/4 1 + m 2 h 1/4 2 = 0 allow both h r 1 and h r 2 to be written in terms of the total unperturbed energy h. Thus, the amplitude can be recast in terms of m 1 , m 2 , ǫ, α m and h as follows: The last identity in (12) is due to the fact that in presence of a single resonance, the motion of the system is tangent to the unperturbed energy surface, and for this particular model such a surface is given by h = h 1 + h 2 , which leads to ∆h 2 = −∆h 1 .
We note that the width of the resonances at O(ǫ, 1/23 2 ) depends on the harmonic numbers in the manner shown in Fig. 1, so that those resonances with values of γ = m 2 /m 1 out of the range [0.5, 2.5] should be narrow.
A glance at both eq. (12) and Fig. 1 reveals that those harmonic vectors with any of its components equal to zero do not change the energies, and consequently, they should not be considered as resonant vectors.
Further, since the perturbation terms are even, if m is a resonant vector then −m is also a resonant one (both corresponding to the same resonance). Hence, for each resonance just one representative resonant vector can be considered, encompassing into its concomitant coefficient the contribution corresponding to its opposite vector as well. All the relevant data required to compute the width of each resonance at O(ǫ, 1/23 2 ) is displayed in Table II, where we have also included the value of h r 1 corresponding to a total unperturbed energy of h = 1/(4β 4 ) ≈ 0.485 (the one used in [10]).
We have computed the resonance widths corresponding to ǫ in the range [0, 0.5]. Fig. 2 displays both the maximum and minimum values of h 1 for each resonance vs. the perturbative parameter ǫ; the total unperturbed energy being h = 1/(4β 4 ). We observe that for ǫ ∼ 0.15, the (6, −1), (4, −1), (2 − 1) and (4, −3) resonances do overlap, but lie far away from the (2, −3) and (2, −5) resonances. Therefrom we could infer that the energy surface presents two   unconnected regions of chaotic motion, so that a global transition to chaos does not take place for ǫ 0.5, which leads to a critical theorical value for the perturbation parameter ǫ c ≫ 0.5.
In those regions of phase space which are far from any primary resonance (i.e. where the diophantic condition holds for every primary resonance), we can introduce new cononical variables, (J, ϕ), in such a fashion that the transformed Hamiltonian consists of a part depending on the new momentum and a perturbation that, though being non-integrable, has an amplitude of O(ǫ 2 ).
where Φ stands for the trigonometric part of the generatrix function of the canonical transformation: After computing the right side of expression (13) one finds: where the coefficients are given by .
There are several relevant facts to be remarked: (i) on working at O(ǫ) we have only considered resonances up to O(1/23 2 ) in the Fourier coefficients; therefore, the series in equation (14) should actually be replaced by a finite sum over those harmonics, m and m ′ , whose associated coefficients (α m and α m ′ , respectively) are such that their product is of order either α 6 1 , α 6 1 /23 or α 6 1 /23 2 ; (ii) there are many different pairs of harmonics (m, m ′ ) at O(ǫ), which combine into the same harmonic n at O(ǫ 2 ); (iii) as a consequence of the evenness of the perturbation term, if n is a resonant vector −n is also a resonant one; (iv) the resonance condition implies that n 1 n 2 ≤ 0; and (v) the condition of being far from resonances at O(ǫ) implies that we must discard all those harmonics which are a multiple of any resonant vector at O(ǫ, 1/23 2 ).
To cope with the situation set up by the issues (ii) and (iii), we have added all the concomitant contributions into a single coefficient D, namely, where the sum extents to all the harmonics (m, is not greater than α 6 1 /23 2 . Taking into account all the above mentioned considerations, the Hamiltonian can be written in the form: Therefore, in the vicinity of a resonant torus J r , and by recourse of the pendulum approximation, we obtain the new resonant Hamiltonian: where Thus, the maximum displacement of the unperturbed action variables is given by (∆J) r ≡ (J − J r ) max = P r n, with P r = 2 (µ|U 0 |) 1/2 , and the widths of the resonances at O(ǫ 2 ) are given by where the plus sign corresponds to n = m + m ′ and the minus sign to n = m − m ′ , m i ought to be written in terms of m ′ i and n i . As a consequence of this last expression, resonant vectors are compelled to have no null components.   The harmonics satisfying all the stated conditions for the perturbation at order ǫ 2 turn out to be just five, which are listed in Table III. They will be referred to as resonant vectors at O(ǫ 2 , 1/23 2 ). We have computed the concomitant resonance widths for ǫ in the range [0, 0.5]. The results are presented in Fig. 3, where also the resonances corresponding to O(ǫ) have been included. Let us recall that the adopted value for the total unperturbed energy is h = 1/4β 4 . From Fig. 3 and Table III we notice that the arising of the (2, −2) resonance connects the two sets of resonances that at O(ǫ, 1/23 2 ) appeared isolated for the considered ǫ range. The remaining resonances at O(ǫ 2 , 1/23 2 ) appear completely overlapped with either set of resonances at O(ǫ, 1/23 2 ). From this plot, we could derive the critical value for the perturbative parameter, ǫ c ≈ 0.28.

VI. NUMERICAL ESTIMATION OF THE CRITICAL VALUE OF THE PERTURBATION PARAMETER
In this section we empirically estimate the value of ǫ c by means of Poincaré Surfaces of Section (SOS). To this aim, we take the intersections on the plane x = 0 (actually |x| < 10 −8 ) whenever p x > 0, for several initial conditions along the y-axis. Fig.4 displays the SOS's corresponding to ǫ = 0.12 (on the left) and to ǫ = 0.14 (on the right), respectively. There we can distinguish the (2, −1) resonance, very close to the last invariant curve that corresponds to the y-axis periodic orbit - (1, 0)   Therefore, it looks like ǫ c lies somewhere in the range (0.12, 0.14). After performing a rather thorough numerical exploration, we have noted that for ǫ = 0.135, several KAM tori do persit, which are shown in Fig.5, where a zoom in the window [0, 0.4] × [0, 0.3] is presented. There, such KAM tori can be clearly distinguished and are seen to definitively separate both chaotic domains in phase space. Nevertheless, this bounded region of chaotic motion, does not involve the resonances we are taking into account to derive the analitical esimation of ǫ c , but high order ones.
Therefore, from experimental means, we may state that ǫ 0.135 is a good lower bound for the critical value of the perturbation parameter.

VII. THE 3D MODEL
Now we will focus on a 3D version of the dynamical system under study. In cartesian coordinates, its Hamiltonian is given by: Let us notice that for a null value of the perturbative parameter we recover the three independent one dimensional quartic oscillators, whose solutions x(t) and y(t) are the ones given by equation (2), while z(t) allows for the expression: where z 0 (h 3 ) = 4βh 3 1/4 . In terms of the action-angle variables of the unperturbed Hamiltonian, the complete Hamiltonian (19) can be recast as: where the new quantitiesV 1j (I) ≡ 2 5/2 3β 4 I 2/3 1 I 1/3 j have been introduced; the ± sign meaning that both terms are included in the series.

VIII. RESONANCES AT O(ǫ)
The perturbation given in eq. (22) shows that for each combination of n, m and k there result 8 harmonic vectors m. Again, due to the even character of the perturbation, we take just one representative resonant vector, m, whose coefficient α m also encompasses the contribution of its opposite vector −m; and keep only those harmonics such that O(α m ) O(1/23 2 ).
Since the harmonic vectors at O(ǫ, 1/23 2 ) can be splitted into two groups, we denote by Y the subset of vectors whose third component is zero and by Z that of vectors having their second component equal to zero. Therefore, the perturbation can be written in the fashion: On applying the resonance condition, m · ω = 0, with m ∈ Z 3 /{0}, to the unperturbed system, the following relation between the actions in each degree of freedom is obtained: so that each resonant vector could not have all its three components of the same sign. Furthermore, whenever a resonant vector has two of its components equal to zero, we get a null amplitude for the perturbation term. Thus, we obtain twelve different resonant vectors, grouped in the following set: For those vectors m ∈ Y, the resonance condition together with the energy conservation condition define a curve in action space (not just a point as in the 2D case) given by: where The concomitant curve for m ∈ Z is given by  . In energy surface both kind of resonances define straight lines. On applying the pendulum approximation, we obtain a similar resonant Hamiltonian to that given by equation (9), namely, where V m =V 12 (I r )α m for m ∈ Y, and V m =V 13 (I r )α m for m ∈ Z.
In energy space, the resonance widths in each degree of freedom, are adequately described by Let us remark that, while in the 2D model the resonance and energy conservation conditions force the resonance width to depend on just one variable, either h, h r 1 or I r 1 , in the 3D model the resonance width is a function of two independent variables, which we have chosen to be I r 1 and h. With the widths computed by means of (28), we can trace the displacements of the resonant energies, h r i + (∆h i ) r m and h r i − (∆h i ) r m , for values of I r 1 ∈ [0, I max ]. Therefore, following [10] we perform the global change of coordinates: where , adopting the value h ≈ 0.485, to finally display in Fig. 6 the region of energy surface occupied by the structure of resonances at O(ǫ, 1/23 2 ) for two different values of the perturbative parameter. Let us remark that many of the resonances in V r (ǫ, 1/23 2 ) are barely observable due to their thinnes and close proximity to a boundary.

IX. RESONANCES AT O(ǫ 2 )
As in the 2D case, we perform a canonical transformation in order to remove the perturbation at O(ǫ), the associated generatrix function being of the form: where The new Hamiltonian is described by the same formal expression given in equation (13), where the first and second terms within braces adopt the values: where δ m,m ′ = 2 3 3 2/3 β 20/3 α m α m ′ and Therefore, the Hamiltonian may be recast as: where A is the set of harmonic vectors arising through any combination of vectors from Y ∪ Z and whose first nonzero component is positive, and the coefficients a, b, c, and d are defined as follows: with J r 1 ∈ [0, ( h A ) 3/4 ]. Let V r (ǫ 2 , 1/23 2 ) be the set of resonant vectors which belong to A and that can be constructed by at least one pair (m, m ′ ) such that O(α m α m ′ ) ≤ O(1/23 2 ). On computing the elements of V r (ǫ 2 , 1/23 2 ) we learn that this set consists of fifteen vectors with one null component together with forty eight having all its three components different from zero. Fig. 7 shows the area of the canonical energy surface (h = 0.485) occupied by those resonances at O(ǫ 2 , 1/23 2 ) that have one null component together with the (2, −1, −1) resonance, for the same two values of the perturbative parameter used for Fig. 6.
Such a subset of V r (ǫ 2 , 1/23 2 ) as well as the complete set V r (ǫ, 1/23 2 ) have been considered in Fig. 8-left for ǫ = 0.005. This picture should be compared with the contour-plot obtained by means of the MEGNO for the very same value of ǫ, which is displayed in Fig. 8-right (taken from [10]). This numerical exploration evinces that the resonances which strongly manifest are those with just one null component (i.e. the straight ones) and the (2, −1, −1) resonance (the one showing a curved shape).
A glance at Fig. 8-left reveals that in some intersections between O(ǫ) and O(ǫ 2 ) resonances, the widths of the latter tend asymptotically to infinity. This is due to the emergence of small denominators in the Fourier coefficients of the perturbation, a fact that reminds us that the canonical transformation perfomed in order to eliminate the perturbation terms proportional to ǫ is no longer valid in the neighbourhood of any O(ǫ) resonance.
A good example of this behaviour is the intersection between the (0, 1, −1) and the (2, −1, 0) resonances. The former is an O(ǫ 2 ) resonance that starts on the top right-hand corner of the energy surface and gets through the middle of it while the latter is of O(ǫ) and can be identified as the widest of the resonances departing from the bottom right-hand corner of the energy surface.

X. ANALYTICAL ESTIMATE OF THE CRITICAL VALUE OF THE PERTURBATION PARAMETER FOR THE 3D SYSTEM
As it can be seen from Fig. 8, many resonances have a triangular shape. Such is the case of all those resonances associated to a vector having either its second or its third component equal to zero.
It can be demonstrated instead that for a resonant vector with its first component n 1 equal to zero, the coefficientes a(n, m, ±m ′ , J) and b(n, m, ±m ′ , J) are null ∀ (m, m ′ ). Thus, the contribution of such a vector proceeds only through its concomitant coefficients c(n, m, ±m ′ , J) and d(n, m, ±m ′ , J). Further, from equation (37) it can be stated that the resonance width tends to zero when J 1 approaches either 0 or (h/A) 3/4 . Consequently, the regions encompassed by the separatrices of resonances for which n 1 = 0 do not have a triangular shape.
On estimating ǫ c for the 3D model, we are compelled to make somewhat strong simplifications: (i) we take as ǫ c the value of ǫ for which the total area covered by resonant regions (A r ) equals 90% of the whole area of the energy surface (A h ); (ii) we approximate by triangles the resonant regions corresponding to resonant vectors with one null component; (iii) we approximate by two triangles the resonant region corresponding to the (2, −1, −1) resonance; (iv) we do not consider any further resonance; (v) we add up the area of each resonance disregarding the intersections due to crossings of resonances, so that those regions corresponding to two different resonances are considered twice.
In Fig. 9-left we have plotted the fraction A r /A h for the perturbation parameter varying in the range ǫ ∈ [0.00001, 0.2]. There it can be observed that A r (ǫ) reaches 90% of A h for some ǫ c between 0.03 and 0.04. This result is in quite good agreement with that arising from Fig. 9-right (taken from [12]) which displays the fraction of chaotic motion according to the MEGNO values, for the same range of ǫ.

XI. DISCUSSION
We have checked out the accuracy of the overlap criterion when applied to a simple near-integrable Hamiltonian system in both its 2D and 3D version. To this end, we have computed the unperturbed resonances up to order O(ǫ 2 ) for both systems, and modelled each resonance by means of the pendulum approximation in order to estimate the theoretical critical value of the perturbation parameter for a global transition to chaos.
By performing several surface of sections for the the 2D case we have derived an empirical value to be compared to our theoretical estimation, both being in good agreement. For the 3D case a theoretical estimate of the critical parameter has been attained, which is shown to match the one given in [12], where such a value is achieved on computing the fraction of chaotic motion vs. ǫ according to the MEGNO values.
Let us remark that the conception of transition to global chaos assumed for the numerical estimate of ǫ c in the 2D case is of a different nature from the one adopted for the 3D system. Actually, the 2D system is considered to be globally chaotic if the chaotic component of phase space appears almost fully connected, while in the 3D case the system is regarded as globally chaotic when at most 10% of the energy surface corresponds to invariant tori. Notice that in the latter case, though it is very likely that the chaotic component be connected when resonances do overlap in a mostly chaotic phase space, nothing could be asserted about the existence of a fully connected region of unstable motion (see [12], [14] and [15] for a thorough discussion).
Therefore, from both theoretical and numerical results we may assert that a suitable estimate for the critical value of the perturbation parameter could be obtained by means of the overlap criterion when considering resonances up to O(ǫ 2 ). Indeed, regarding terms just up to O(ǫ) largely overestimates ǫ c , as already shown by [1] for the Standard Map.