Systematic search for islets of stability in the standard map for large parameter values

In the seminal paper (Phys. Rep. 52, 263, 1979), Boris Chirikov showed that the standard map does not exhibit a boundary to chaos, but rather that there are small islands (islets) of stability for arbitrarily large values of the nonlinear perturbation. In this context, he established that the area of the islets in the phase space and the range of parameter values where they exist should decay following power laws with exponents -2 and -1, respectively. In this paper, we carry out a systematic numerical search for islets of stability and we show that the power laws predicted by Chirikov hold. Furthermore, we use high-resolution 3D islets to reveal that the islets volume decays following a similar power law with exponent -3.


Introduction
The standard map, also known as Chirikov-Taylor standard map, is a paradigmatic 2D area-preserving map given by [1] where p n and θ n are taken modulo 2π; and K > 0 is a constant whose role is to increase the nonlinear perturbation.The physical meaning of θ and p depends on the system under consideration.This is a quite general model which describes the dynamics of a nonlinear oscillator under periodic perturbations, and it has been broadly studied in the context of nonlinear dynamics.Some real-world physical situations described by the standard map include a kicked rotor [2] and the motion of a charged particle in a magnetic bottle [1].
For small values of K (say K < 4), the phase space (θ, p) is occupied by a chaotic sea and a main Kolmogorov-Arnold-Moser (KAM) island that covers a significant area of the phase space.It is noteworthy that all the KAM islands in this system surround stable periodic orbits.The coexistence of stable and chaotic dynamics is a hallmark of Hamiltonian systems and, in words of Zaslavsky, "is one of the most striking and wonderful discoveries ever made" [2].As K increases, the area of the KAM island is irregularly reduced until being eventually destroyed for K ≈ 7.For higher values of the parameter K, a fully chaotic regime might be expected.However, small regions ("islets") of stability emerge in an approximate periodic manner, even for arbitrarily large parameter values [3].Akin to KAM islands, islets of stability are regions filled with KAM tori where orbits are dynamically trapped.However, islets are not directly related to the main KAM island.As a matter of fact, they do not emerge in bifurcations of the main family of periodic orbits.Instead, they appear in saddle-node bifurcations out in the chaotic sea [4].Although islets have not received much attention in the literature, their existence is not a unique feature of the standard map, but they also appear in general Hamiltonian systems and area-preserving maps [5,6,7].
Chirikov himself suggested that the length of the parameter interval ∆K where the islets exist and their area A in the phase space should decay following power laws like ∆K ∝ K a1 and A ∝ K a2 .The predicted exponents are a 1 = −1 and a 2 = −2.The latter exponent was numerically investigated by Contopoulos et al. in 2005 [8].They classified the islets in two different types and obtained a decay exponent (−2.45 and −2.3) for each type.Although these numbers are relatively close to the theoretical value −2 (relative error about 15 − 20%), they do not provide definitive numerical evidence supporting the idea that these laws govern the decay of magnitudes in the system.To the best of our knowledge, there have not been more refined estimations of these exponents.
In this paper, we conduct rigorous numerical simulations to validate the accuracy of Chirikov's predictions, which we find to be remarkably precise.In particular, we obtain the predicted decay exponents with relative errors of less than 0.6%.Moreover, we extend the previous predictions to the case of the volume in the (θ, p, K) space, obtaining an exponent −3.
The manuscript is organized as follows.First, in Section 2, we introduce a systematic approach for computing islets of stability.In Section 3, we explore numerous islets across a wide range of parameters, demonstrating that their length, area, and volume strictly follow power laws.Finally, in Section 4, we summarize the key findings and discuss the potential implications of the results.

Islets of stability
A general overview of the standard map dynamics can be obtained simply by plotting orbits in the phase space.However, since we aim to locate islets of stability in a systematic manner, we need a method to detect the last KAM curve, i.e., the boundary between a KAM island and the chaotic sea.In this section, we present the methods we have implemented to accurately locate the islets.
Due to the modulo operation, the standard map is a mapping from the two-dimensional torus to itself.Therefore, any orbit is confined to the phase space region [0, 2π] × [0, 2π].Nevertheless, orbits within the chaotic sea will eventually visit the entire phase space, except the region bounded by the last KAM curve.Using this simple fact, and following Sanjuán et al. [9], we define regions ("leaks") in the chaotic sea.If an orbit never falls into them, we classify it as a part of a KAM island.Naturally, for this method to work the leaks must be located in an area that is known a priori to be in the chaotic sea.In particular, we define a left leak L 1 ≡ [0.1π, 0.3π] × [0, 2π] and a right leak . This choice introduces two rectangular leaks of width 0.2π that are symmetrical about θ = π.We highlight that this choice is arbitrary and any region in the chaotic sea can be used to serve this purpose.
Once the leaks are defined, we divide the phase space into a grid of initial conditions and we iterate the system for each of them.If after one or more iterations an orbit falls into one of the leaks, we stop the iteration process and we classify the corresponding initial condition as a part of the chaotic sea.On the contrary, if after a fixed maximum number of iterations the orbit has never fallen into the leaks, we classify the corresponding initial condition as a part of a KAM island.To visualize the result, we assign blue (red) color to the initial conditions whose orbits fall into the left (right) leak.Finally, we assign white color to initial conditions belonging to a KAM island.
The procedure that we have just introduced is equivalent to the numerical simulation of exit basins in open Hamiltonian systems [10].For this reason, we will refer to the corresponding data and plots as "exit basin diagrams".Some examples computed for different values of K are depicted in Fig. 1.For K < 4, the main KAM island surrounds a periodic orbit located at (θ, p) = (π, 0) [or, equivalently, (θ, p) = (π, 2π)], as shown in panels (a) and (b).For higher values of K [see panels (c) and (d)], the periodic orbit moves along two symmetry lines p = 2θ and p = 2θ−2π, which are represented with gray dashed lines in Fig. 1.Some periodic orbits generating islets of stability also appear over these symmetry lines.
Since islets appear in a reduced range of parameter values, it is difficult to identify them by computing exit basin diagrams in the phase space.Furthermore, as we shall show later, the interval ∆K where an islet exists is reduced as K increases.Thus, this task becomes even more arduous for large parameter values.However, the fact that the main stable periodic orbits, and hence the KAM islands, lie on symmetry lines can help us to narrow our search.In particular, we can preliminarily locate islets by constructing exit basin diagrams in which we choose initial conditions along a given symmetry line for different values of K.As an example, in Fig. 2 we show the exit basin diagram computed along the symmetry lines p = 2θ − 2π and p = 2θ for parameter values in the range K ∈ [0, 20].We emphasize that other symmetry lines, such as θ = p, could be used to find different islets.For low values of K, we observe the main KAM island, which is destroyed for K ≈ 7.This island possesses a characteristic fractal tree-like structure that is ubiquitous in area-preserving maps and Hamiltonian systems (see, for example, Refs.[6,11]).After the destruction of the main island, the fully chaotic regime is reached, but only momentarily.The first islet of stability appears for K ≈ 9.19, and it is followed by a second islet appearing for K ≈ 11.85.As it has been shown by Contopoulos et al. [8], to precisely study the characteristics of the islets it is necessary to classify them in different types.Each islet type is characterized by the period of its central periodic orbit and can be found in a particular symmetry line.In this work, we focus our attention on islets surrounding period-4 orbits (type I) and period-2 orbits (type II).These periodic orbits, which are periodic on the torus but nonperiodic on the cylinder, have been referred to as "accelerator modes" in the literature [12], and their diffusion properties have received some attention [13,14,15].The first islet of each type, which are also highlighted with white dashed rectangles in Fig. 2, are shown with high resolution in panels (a) and (c) of Fig. 3.As can be noticed by looking at the scale of the figures, type II islets are rather bigger than type I.
As K is increased, a scaled copy of each type of islet appears in approximate parameter intervals of 2π.Their size is reduced as K increases, but their shape remains unaltered.In other words, exact copies of the islets appear recurrently, simply changing their size as the leaves in a fern.This can be intuitively seen in Fig. 2, where it is clear that the second islets of each type are smaller.To illustrate the astonishing similarity between islets of the same type, we depict islets of types I and II for values of K close to 1000 in Fig.  and Fig. 3d.We can see with naked eye that their main features, ignoring their size, correspond to a great extent to the islets appearing for low values of K (panels (a) and (c) of Fig. 3).

Decay of magnitudes
In this section, we compute the length, area and volume of the islets, as well as we obtain the exponents governing the decay of these magnitudes.To accomplish this, we systematically locate multiple islets following the procedure described in Sec.II.Afterward, we construct high-resolution 3D exit basin diagrams in the (θ, p, K) space.As an illustrative example, Fig. 4 displays the first islet of each type.
These 3D representations offer a comprehensive view of the complex evolution of the islets as the parameter K is increased.The rich fractal structures of the islets can be seen with great detail by zooming in the figure.From these data, a Poincaré surface of section for a fixed value of K brings exit basin diagrams equivalent to those showed in Fig. 1.In a similar fashion, defining the Poincaré section p = 2θ − 2π, the exit basin diagrams of Fig. 3 are obtained.
The length, area and volume of an islet can be obtained from a 3D exit basin diagram.The main ingredients of the procedure we have followed are sketched in Fig. 5.In this figure, the gray volume corresponds to the same type II islet plotted in Fig. 4b.On the other hand, we highlight specific cross sections for K = const (indicated by the blue areas within the 3D islet), which represent islets in the phase space (θ, p).The area A of an islet in the phase space evolves as K increases, as shown on the right blue curve.Therefore, to characterize different islets in terms of their area, it is necessary to establish a reference point.In our case, we have chosen as a reference their maximum area, max(A).Finally, the length ∆K of an islet is simply the parameter range where it exists.
According to Chirikov, the length and maximum area of the islets should follow power laws of the form where a 1 = −1 and a 2 = −2.From the perspective of numerical simulations, K is the parameter value where an islet appears (see K i in Fig. 5).
As far as we know, no power laws for the volume have been suggested in the literature.However, given that the length decays as 1/K and the area as 1/K 2 , it leads us to conjecture that the volume should decay as 1/K 3 .Therefore, a similar power law with exponent a 3 = −3 is expected.To show that these decay laws hold, we have studied all the islets appearing in the range of parameter values K ∈ [9, 500], where we found a total of 157 islets (79 of type I and 78 of type II).For each islet, we have computed a 3D exit basin diagram with a 1000×1000× 1000 resolution.
Once the length, maximum area, and volume have been calculated, we estimate the exponents a 1 , a 2 , and a 3 by fitting least squares lines on log − log plots.These logarithmic plots, along with the corresponding least squares lines, are depicted in Fig. 6.Detailed results, including the decay exponents a, linear correlation coefficients r, and relative errors δ, are presented in Table 1.Table 1 Slope a and linear correlation coefficient r for the straight lines depicted in Fig. 6.The slope of the lines corresponds to the exponent of the power laws given by Eq. ( 2) and Eq. ( 3).The relative error δ is calculated as |(ā − a)/ā|, being ā the theoretical value and a our estimation.
Fig. 6 Logarithmic plots for (a) length, (b) maximum area, and (c) volume of the islets.Blue dots represent type I islet data, while red dots represent type II islet data.Following the same color scheme, each data set is accompanied by a straight line obtained using the least squares method.The slope of the lines corresponds to the exponent of the power laws given by Eq. ( 2) and Eq. ( 3).
The experimental values match well with the theoretical predictions.Notably, the relative errors range from 0.0017% to 0.57%.Additionally, the linear correlation coefficients are very close to −1 in all cases, suggesting a robust linear relationship between the variables.

Conclusions and discussion
In summary, in this paper we leverage the symmetries of the standard map to systematically identify islets of sta-bility for a wide range of parameter values.We have accurately determined the exponents governing the power laws describing the decay of the length, maximum area, and volume of the islets.Theoretical values for the first two exponents were already established in the literature by Chirikov, but a rigorous numerical proof had not yet been provided.
In addition to these results, we have carefully described the numerical schemes used to locate islets of stability.These methods can serve as useful tools for future research endeavors aimed at identifying and studying islets in different systems.
Despite their small size, islets of stability can exert a significant influence in a system.As a matter of fact, even small KAM islands can generate anomalous diffusion [16,17] and modify global properties such as decay correlations [18].In certain systems, maintaining stable motion may be preferable over chaotic behavior.In this situation, islets constitute a valuable subset even for large parameter values where the system is mostly (but not fully) chaotic.As a result, to find islets can be advantageous for control schemes or to define strategic parameter values.One particularly interesting potential application arises in the context of plasma confinement in tokamaks [19].As an inkling of such potential application, we have conducted preliminary simulations on the Boozer map [20,21], which is a model for the configuration of the magnetic field lines of a tokamak, and we have verified that it exhibits islets of stability.

Fig. 1
Fig. 1 Exit basin diagram in the phase space showing KAM islands (white color) for parameter values (a) K = 2, (b) K = 4, (c) K = 5, and (d) K = 6.Blue (red) color represents chaotic orbits that fall into the left (right) leak.The gray dashed lines correspond to the symmetry lines p = 2θ and p = 2θ − 2π.

Fig. 2
Fig. 2 Exit basin diagram showing in white color the main KAM island and some islets in the (θ, K) plane.Blue (red) color represents chaotic orbits that fall into the left (right) leak.To obtain these data, once an initial condition θ is chosen, the remaining coordinate is given by p = 2θ − 2π if θ > π or by p = 2θ if θ ≤ π.The islets highlighted with white dashed rectangles are shown with high resolution in Fig. 3 [panels (a) and (c)]. 3b

Fig. 3
Fig. 3 Exit basin diagram in the (θ, K) plane showing islets of stability (white color) of (a,b) type I and (c,d) type II.It is clear that islets of the same type are identical scaled copies.Blue (red) color represents chaotic orbits that fall into the left (right) leak.To obtain these data, once an initial condition θ is chosen, the remaining coordinate is given by p = 2θ − 2π.

Fig. 4
Fig. 4 Islets of stability of (a) type I and (b) type II in the (θ, p, K) space.White color represents initial conditions within an islet, while blue (red) color is used on three of the cube's faces to represent chaotic orbits that fall into the left (right) leak.To properly visualize the islets, the remaining chaotic initial conditions are left transparent.Great detail of the islets arises when zooming in.

Fig. 5
Fig.5Main ingredients of our method for determining the length, maximum area, and volume of the islets.The gray volume corresponds to a 3D islet, while blue areas (cross sections at constant K) illustrate the islet's area for various values of K.The blue curve on the right displays the evolution of the area in the phase space (θ, p) as a function of K. Some characteristic values of the area, including its maximum, are indicated with arrows.The length of the islet is ∆K = K f − K i .