Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://chaos.phys.msu.ru/loskutov/PDF/Reprint_JPA_2011_44_175102.pdf Äàòà èçìåíåíèÿ: Fri Apr 8 14:08:53 2011 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:19:11 2012 Êîäèðîâêà:

IOP PUBLISHING J. Phys. A: Math. Theor. 44 (2011) 175102 (12pp)

JOURNA

L OF

PHYSICS A: MAT

HEMATICAL AND

THEORETICAL

doi:10.1088/1751-8113/44/17/175102

A family of stadium-like billiards with parabolic boundaries under scaling analysis
´ Andre L P Livorati1 , Alexander Loskutov1 ,2 and Edson D Leonel1
1

´ Departamento de Estat´ istica, Matematica Aplicada e Computacao, Univ. Estadual Paulista, ¸~ Av.24A, 1515, Bela Vista, CEP: 13506-900, Rio Claro, SP, Brazil 2 Physics Faculty, Moscow State University, Moscow 119992, Russia E-mail: andrelivorati@gmail.com

Received 8 December 2010, in final form 3 March 2011 Published 4 April 2011 Online at stacks.iop.org/JPhysA/44/175102 Abstract Some chaotic properties of a family of stadium-like billiards with parabolic focusing components, which is described by a two-dimensional nonlinear areapreserving map, are studied. Critical values of billiard geometric parameters corresponding to a sudden change of the maximal Lyapunov exponent are found. It is shown that the maximal Lyapunov exponent obtained for chaotic orbits of this family is scaling invariant with respect to the control parameters describing the geometry of the billiard. We also show that this behavior is observed for a generic one-parameter family of mapping with the nonlinearity given by a tangent function. PACS numbers: 05.45.-a, 05.45.Pq, 05.45.Tp (Some figures in this article are in colour only in the electronic version)

1. Introduction A two-dimensional billiard consists of a domain Q with a piecewise-smooth boundary Q inside which a point mass (billiard ball) moves freely along straight lines with a constant velocity. Reaching the boundary, the particle is reflected from it elastically such that the normal component changes its sign while the tangent one is preserved. This condition leads to the following consequence: the incidence angle is equal to the angle of reflection [1]. In the contemporary form, the notion of billiard is known since Birkhoff [2] who issued the challenge of the billiard ball motion a manifold with an edge. A more complete analysis regarding the mixing property in hard ball systems has been performed by Krylov [3]. However, thanks to the papers of Sinai [4, 5], who developed some mathematical ideas for the study of ergodic properties for a dispersing planar billiards, a new epoch in the rigorous analysis of billiards began.
1751-8113/11/175102+12$33.00 © 2011 IOP Publishing Ltd Printed in the UK & the USA 1


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

(a)

(b)

Figure 1. A stadium-like billiard and its geometric parameters. (a) Bunimovich stadium with fully chaotic properties. (b) A billiard corresponding to a nearly integrable system (b l ).

These ideas led to an avalanche-like growth of physical and mathematical investigations of billiard models. A remarkable advance has been achieved in a solution, in some form concerning the Boltzmann ergodic hypothesis for the dynamics of hard ball gases [4]. Nowadays the study of billiards relates to such fields as non-equilibrium statistical mechanics [1], hyperbolic systems [6, 7], optics [8, 9] and several physical experiments such as microwave billiards [10, 11], ultra-cold atoms trapped in a laser potential [12­15] and also mesoscopic quantum dots [16]. All planar billiards can be subdivided into dispersing, focusing and neutral. Dispersing, sometimes also called as Sinai billiards, consists of both dispersion and neutral components. One popular example of such billiards is Lorentz gas [1, 4]. On the basis of the analysis of a 2D Lorentz gas, it has been proved that the dynamics of purely deterministic systems may be ergodic with mixing and similar to the Brownian motion [17]. The boundary of focusing billiards consists of only focusing components or focusing components which are connected by straight lines. Typical examples of such billiards belong to Bunimovich [18, 19]. This billiard includes two arcs and two parallel segments joining them (see figure 1) and therefore looks like a stadium. For such a billiard it was proved that for some geometric parameters a, b and l it possesses the property of ergodicity and mixing. Conditions of the chaoticity in 2D plane billiards are described in [1]. The chaos generated in dispersing and focusing billiards is different. In the former case, two narrow beams of billiard particles being reflected from a dispersing component will be diverged. In the latter case, after the reflection, two parallel beams converge in a focusing point. Chaos in such billiards appears if the time, during which the beam is diverged, is larger than the convergence time. Thus, depending on the billiard geometry one can obtain a principal difference in the particle dynamics. In general, dynamics of the billiard particle may be of three qualitatively different kinds: (i) completely integrable, i.e. regular; (ii) ergodic with mixing; and (iii) depending on initial conditions, regular or mixing. The regular dynamics is inherent in billiards with the boundary in the form of a circle and ellipse [1]. Sinai billiard and Bunimovich stadium, see figure 1(a), possess ergodicity and mixing properties [1, 5, 17, 19]. The third kind can be observed in billiards of many configurations (see, e.g., [7, 20, 21]). In particular, the Bunimovich stadium at certain parameter values (figure 1(b)) has both regular and mixing particle dynamics. For such billiard systems, a set of Komolgorov­Arnold­Moser (KAM) islands can be observed. They are usually surrounded by chaotic layers and invariant spanning curves (also called invariant tori), separating different regions of the phase space.
2


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

Often, the dynamics of the system depends on the control parameters. Particularly, such parameters can control the intensity of the nonlinearity of the system and, as they change, average quantities of the dynamics can experience subtle modifications, therefore characterizing a phase transition [22]. Many of these transitions can be described using scaling arguments where critical exponents play an important role in the description of the average properties. In particular, a phase transition from integrability to non-integrability can be studied [8], and therefore classes of universality can be obtained [23]. The scaling formalism has been shown valid to describe several different kinds of phase transitions in billiard systems considering both static [8], time-dependent non-dissipative [24, 25] as well as dissipative cases [26­28]. In this paper, we consider a family of stadium-type billiards with fixed parabolic boundaries. We show that geometric parameters of this family can be rescaled in the sense that the maximal Lyapunov exponent of chaotic components of the billiard map is scaling invariant under the billiard geometry. The organization of the paper is given as follows. In section 2, we carefully discuss the construction of the map which describes the dynamics of billiard particles. Description of a transition from local to fully chaotic dynamics is also presented and is based on Bunimovich's defocusing mechanism. Section 3 presents our numerical results related to calculations of the maximal Lyapunov exponent. A discussion concerning scaling invariance is given here as well as some numerical results for the deviation of the average angle between the trajectory and the vertical line at the collision point. We also give a generalization of the mapping and discuss similar transition for a one-parameter family of mapping. Final remarks and conclusions are drawn in section 4.

2. A family of maps for stadium-like billiards In this section, we discuss in detail how to construct the corresponding family of maps that describe the dynamics of the billiard particle. Suppose that the stadium family has the shape shown in figure 1(b). Then, for b a l , the billiard loses the K-property and turns out to be a nearly integrable system. In this case, its phase space consists of stable fixed points surrounded by KAM islands where the dynamics of the particle is regular, and a chaotic sea that is characterized by a positive Lyapunov exponent. As usually considered in the literature, to analyze the billiard system, we take into account only collisions with the focusing components, i.e. consider the unfolding of the billiard table (see figure 2). Let us introduce the coordinates as shown in this figure. The motion of the billiard particle inside the billiard area Q (see figure 1) is given by a map (xn+1 ,n+1 ) = T(xn ,n ), where xn is a collision point taken as mod a , and n is the angle between the particle trajectory and the vertical line. The angle is measured counterclockwise. Suppose now that the focusing components of the billiard boundary are parabolic. Moreover, assume that they are described by a quadratic function of a general type f(x ) = Ax 2 + Bx + C, (1)

where A, B and C are constants and can be determined by boundary conditions. Considering the case where f(0) = 0 and f(a ) = 0, it is easy to see that C = 0 and B = -Aa . For the case of f(a /2) = -b, one can obtain that A = 4b/a 2 and B = -4b/a . Replacing the
3


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

Figure 2. Coordinates in a stadium unfolding and its dynamical variables.

expressions obtained from the boundary conditions above in equation (1), we arrive at the equation for the parabolic boundary: (x ) = 4bx (x - a), a2 (2)

which describes the equation of the lower part of the boundary. Applying the procedure in a similar way, but taking into account the negative sign in the quadratic term and a vertical translation of l, one can easily obtain the equation of the upper boundary. Since we know the analytical expressions of the boundary, the next step is to obtain the map T. From geometrical considerations, as shown in figure 2, and assuming that l b, we can see that xn+1 = xn + l tan n . (3)

We would like to emphasize that the condition l b, corresponds to an approximation of the real billiard model. To obtain n+1 , we first define an auxiliary function (x ) = arctan (x ), which is the slope of a tangent line at the collision point with the boundary. Thus, expanding (x ) in Taylor series, taking the lower order term and taking into account equation (2), we have (x ) 4b (2x - a). a2

Considering that the angles and are obtained from the normal line as shown in figure 2, one can see that n+1 + + n+1 = , 2 and + = + = /2. From these expressions we found that n+1 + n+1 = . Therefore, n+1 = 2 - n . Defining a dimensionless variable = x/a with [0, 1) and using equation (3), we obtain the known map for the stadium-like billiard with static parabolic
4


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

(a)

(b)

(c)

(d )

Figure 3. Time evolution of a single initial condition for the static stadium-like billiard with parabolic boundaries. The control parameters used were a = 0.5, l = 1 and (a) b = 0.07; (b) b = 0.05; (c) b = 0.01 and (d) b = 0.005. In particular for (c) we have the fixed points of period 1 labeled as circles, period 2 represented as squares, and period 3 identified as diamonds.

boundary, that is the same map obtained in [29]: = + l tan mod 1 n+1 n n a T: = - 8b (2 - 1). n+1 n n+1 a

(4)

Since the stadium billiard is axial-symmetric, the transformation T(n ,n ) = (n+1 ,n+1 ) should have a symmetry with respect to 1 - and - . Therefore, it is enough to obtain results for the region of a non-negative value of . Figure 3 shows the evolution of one initial condition chosen in the chaotic sea, for very long time in the phase plane for different control parameters. In particular for figure 3(c), we show evolution of some initial conditions chosen near the fixed points. The same procedure shown in figure 3(c) can be extended to figures 3(b) and (d). One can see that as the parameter b assumes large values, the chaotic component of the phase space is larger. However, as the control parameter b is reduced, keeping fixed both l and a, more stable islands appear in the phase space. Also, invariant curves appear corresponding to regular oscillations. Thus, a phase transition from global chaos to partially chaotic dynamics is observed.
5


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

As is known [30], the transition to fully chaotic dynamics takes place when all the stable fixed points lose their stability; thus, only unstable orbits are present in the phase space. Therefore, in order to apply the stability criteria we must first obtain the fixed points for the map given by equation (4). Their expressions are n = 1/2+ m and n = arctan(ma / l ), which are indeed elliptical points. The parameter m Z corresponds to the number of mirrored stadiums that the particle went through (using the unfolding method) before colliding with the other focusing boundary. The stability of the period-1 fixed points is obtained by the eigenvalues of the Jacobian matrix: n+1 n+1 n n , J =
n+1 n n+1

where

n n+1 l = , n a cos2 n

n+1 = 1, n

n+1 n+1 16b 16bl , =- =1- 2 . n a n a cos2 n Since det J = 1, the map preserves the area in the phase space. If | Tr J | 2, then all the stable fixed points bifurcate and become unstable. This leads to 4bl > 1. (5) a2 If the inequality given in (5) holds, then global chaos, in the sense that no stable periodic fixed points exist, is observed in the system described in equation (4). In billiards which have focusing components of a constant curvature, the particle dynamics is chaotic if all these components complemented up to a circle belong to the billiard table Q [31]. The inequality given in (5) is approximately the same condition where the Bunimovich defocusing mechanism [1, 31] takes place for the classical stadium billiard (see figure 1(a)). This mechanism basically consists of a narrow parallel beam of rays, that after being focused due a reflection with the focusing boundary may pass through a focusing point and become divergent provided that a free path between two consecutive reflections from the boundary is long enough. The mechanism of defocusing works under the condition that divergence prevails over convergence. 3. Numerical results We start the section by discussing our numerical results for the positive Lyapunov exponent of chaotic components of the phase space. The Lyapunov exponent has been widely used to quantify the average expansion or contraction rate for a small volume of initial conditions. If the Lyapunov exponent is positive, the orbit is said to be chaotic. On the other hand, a non-positive Lyapunov exponent indicates regularity and the dynamics can be in principle periodic or quasi-periodic. The Lyapunov exponents are defined as follows [32] (see for example [33, 34] for applications in higher dimensional systems): 1 j = lim ln n , j = 1, 2, j n n n where j , are the eigenvalues of the matrix M = n=1 J(i ,i ) and Ji is the Jacobian matrix i evaluated over the orbit. Figure 4 shows the behavior of the average value of the positive
6


J. Phys. A: Math. Theor. 44 (2011) 175102
6
2.8

A L P Livorati et al

5 4 3

2.4 2 1.6 0.06 0.062 0.064 0.066

Numerical Data



2 1 0 -3 10

bc

bc
10
-2

10

-1

10

0

b
¯ Figure 4. Plot of â b. The sudden change characterizes the transition from chaos to stability. The control parameters used were a = 0.5 and l = 1. The zoom-in shows the transition better, where b = 10-4 .

Lyapunov exponent as a function of the control parameter b of the map given in equation (4). For this plot, we kept fixed l = 1 and a = 0.5. Each point in the figure was obtained for seven different initial conditions randomly chosen along the chaotic component iterated up to 108 collisions with the boundary (i.e. iterations for the map in equation (4)). We emphasize that no more statistics are needed because of the excellent convergence of the positive Lyapunov exponent as a function of time. To comment on this, given an initial condition in the chaotic sea, the time evolution of the Lyapunov exponent fluctuates for short time and then converges to a constant plateau for about 104 iterations of the mapping. It also stays constant for larger iterations. Other initial conditions chosen in the chaotic sea exhibit similar behavior therefore converging then to a very near plateau of the previous initial condition within a small uncertain numerical range. Also, long time simulations of collisions assure a good convergence of the algorithm. The error bars in figure 4 correspond to the standard deviation of the seven samples. Similar procedures analyzing transitions as a function of the maximal Lyapunov exponent were considered previously [35, 36]. ¯ We see that there is a sudden change in when the transition to global chaos occurs (absence of periodic and stable fixed points), i.e. when the defocusing mechanism starts to work. For fixed a and l, the critical value, say bc, where this transition occurs is bc = 0.0625. The zoom-in shown in figure 4 corresponds to the amplification of the ¯ behavior of . The step used in the parameter b was b = 10-4 , leading to a very accurate investigation. ¯ The sudden jump in the behavior of is also observed for other different combinations of the control parameters. If the control parameter l is raised with a fixed, according to the condition given in equation (5), the transition is observed for a smaller b. Figure 5(a) ¯ shows several curves of â b for different values of l. All curves exhibit qualitatively ¯ the same general behavior; however, the sudden transition on occurs at different critical parameters bc. For the limit b 0, the Lyapunov exponent tends to zero, thus leading to the regular case. As far as we have observed, the transition is smooth and the Lyapunov exponent tends to zero monotonically.
7


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

10

0



10
(a)

-1

l l l l
-6

=1 =10 =100 =1000

10

10

-4

10

-2

10

0

b
0

10 10

-1

10

-2

a = 0.1 a=1 a = 10


10
(b)
-3

10

-4

10

-2

10

0

10

2

b

¯ Figure 5. Plots of â b. The control parameters used were (a) a = 0.5and l = 1, l = 10, l = 100 and l = 1000; (b) l = 1 and a = 0.1, a = 1and a = 10.

One can ask about the transition as a function of the billiard parameter a. Indeed, the ¯ same general behavior is observed. Figure 5(b) shows a plot of â b for a fixed l = 1 and three different values of a, as labeled in the figure. The sudden change is evident. ¯ The behavior shown in both figures 5(a) and (b) indicates that seems to be a scaling invariant. In fairness, after rescaling the axis properly by bl /a 2 , a universal representation for ¯ all curves is obtained as shown in figure 6. We see that all curves of , obtained for different control parameters, merge together in a perfect collapse of a single and universal plot. This ¯ is a clear confirmation that is a scaling invariant with respect to the transformation given in equation (5). We can also characterize the behavior of the deviation of the average angle, called , near the critical region where the defocusing mechanism holds. We define as (n, l , a , b) = 1 M
M

i2 (n,l ,a ,b) - i (n,l ,a ,b),
i =1

2

(6)

where M corresponds to an ensemble of different initial conditions while the average angle is obtained by (n, l , a , b) = 1 n
n

i .
i =1

(7)

Equation (6) was iterated up to 108 collisions of the particle with the boundary of the billiard, using an ensemble of M = 5 â 103 different initial conditions chosen along the
8


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

10

0


10
-1

10

-2

10

-3

10

-2

10

-1

bl/a

2

10

0

10

1

¯ Figure 6. Universal plot for different curves of as a function of the rescaling parameter bl /a 2 .

Figure 7. How the defocusing mechanism affects the behavior of curves, characterizing the transition from fully to partially chaotic dynamics. The control parameter b is labeled in the figure while a = 0.5and l = 1.

chaotic sea. Figure 7 shows the behavior of as a function of n for different values of the control parameter b, as labeled in the figure. Since the geometrical parameters are all scaling invariant, we decided to vary only the parameter b, thus keeping a and l constant. We see in figure 7 three different kinds of behavior according to the control parameter. When b is chosen such that 4bl /a 2 > 1 (absence of stable periodic fixed points), has an initial growth and, after some a short number of collisions, bends toward a regime of saturation for large values of n. For such behavior, is restricted to a range (0.8, 1). For this situation, the defocusing mechanism is working.
9


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

¯ Figure 8. Behavior of â K . A zoom-in shows amplification of the sudden jump, near the criticality characterizing the transition from fully to partially chaotic dynamics as a function of the rescaling parameter K.

When b is set to give 4bl /a 2 1, the curves of start to grow for the regime of small n = and then have a sudden decay, thus converging to a constant plateau in the range (0.1, 0.3). Finally, when 4bl /a 2 < 1 (defocusing mechanism is not working), the curves are almost constant for all values of n. Moreover, they are all confined in the range below 0.1 implying that the average angle is very small. Let us now discuss a little more deeply a possible reorganization of mapping (4) in a generic family of one-parameter mapping. Defining conveniently the parameter pn+1 = 2n+1 - 1 and In+1 = (8b/a )pn+1 and setting K = 16bl /a 2 , the mapping (4) is rewritten as T: In+1 = In + K tan(pn ) . pn+1 = pn - In+1 (8)

For this new mapping, all period-1 fixed points become unstable for K > 4 [37]. A plot ¯ of â K is shown in figure 8. A sudden jump is also evident for K 4 which confirms = the sudden increase in when the stable fixed points lose their stability, thus leading the dynamics to be fully chaotic. The zoom-in shown in figure 8 shows such a transition better. In principle, it is expected that such a generalization can also be applied for a whole family of one-parameter mappings with a nonlinearity given by a tangent function.

4. Concluding remarks In general, the universality phenomenon is inherent in wide kinds of dynamical systems where chaos is observed. As is known, a quite natural idea lies at the heart of the universality: certain properties of dynamical systems are completely repeated at several scalings. Specifically, the scaling procedure allows us to understand the integrability­chaos transition. Therewith, the staggering aspect of this approach is quantitative descriptions. The scaling formalism may also provide a useful implement in the characterization of billiard systems. In particular, this technique is useful for the explanation of their chaotic properties. It makes it possible to generalize and describe a quite large class of the billiard models for which chaos is inherent. In this paper, we have investigated the behavior of the maximal Lyapunov exponent for chaotic orbits along the chaotic sea in a stadium-like billiard with parabolic boundary.
10


J. Phys. A: Math. Theor. 44 (2011) 175102

A L P Livorati et al

A transition from stability to chaos was discussed based on the Bunimovich defocusing mechanism. A sudden transition was observed for the positive Lyapunov exponent near the criticality where the defocusing mechanism works. The behavior of the Lyapunov exponent was shown to be scaling invariant with respect to the control parameters of a family of stadiumlike billiards with parabolic components. We demonstrated that the increase of chaoticity in such a family obeys certain scaling rules. In other words, under rescaling of the billiard parameters a universal behavior of the maximal Lyapunov exponent is observed. This sudden jump was also observed for a generic rewriting of the stadium mapping leading to a family of one-parameter map.

Acknowledgments The authors acknowledge fruitful discussions with Professor Vanderlei Marcos do Nascimento. ALPL and AL thank FAPESP. EDL acknowledges support from FAPESP, CNPq and FUNDUNESP, Brazilian agencies. AL also expresses thanks to DEMAC for the kind hospitality during his stay in Brazil.

References
[1] Chernov N and Markarian R 2006 Chaotic Billiards vol 127 (Providence, RI: American Mathematical Society) [2] Birkhoff G D 1927 Dynamical Systems (Am. Math. Soc. Colloquium Publ. vol 9) (Providence, RI: American Mathematical Society) [3] Krylov N 1979 Works on the Foundations of Statistical Physics (Princeton, NJ: Princeton University Press) [4] Sinai Y G 1963 Sov. Math. Dokl. 4 1818 [5] Sinai Y G 1970 Russ. Math. Surv. 25 137 [6] Robnik M and Berry M V 1985 J. Phys. A: Math. Gen. 18 1361 [7] Berry M V 1981 Eur. J. Phys. 2 91 [8] Leonel E D 2007 Phys. Rev. Lett. 98 114102 [9] Persson E et al 2000 Phys. Rev. Lett. 85 2478 [10] Stein J and Stokmann H J 1992 Phys. Rev. Lett. 68 2867 [11] Stokmann H J 1999 Quantum Chaos: An Introduction (Cambridge: Cambridge University Press) [12] Milner V et al 2001 Phys. Rev. Lett. 86 1514 [13] Friedman N et al 2001 Phys. Rev. Lett. 86 1518 [14] Andersen M F et al 2004 Phys. Rev. A 69 63413 [15] Andersen M F et al 2006 Phys. Rev. Lett. 97 104102 [16] Marcus C M et al 1992 Phys. Rev. Lett. 69 506 [17] Bunimovich L A and Sinai Y G 1981 Commun. Math. Phys. 78 479 [18] Bunimovich L A 1974 Funct. Anal. Appl. 8 254 [19] Bunimovich L A 1979 Commun. Math. Phys. 65 295 [20] Oliveira D M F and Leonel E D 2010 Commun. Nonlinear Sci. Numer. Simul. 15 1092 [21] Koiller J et al 1996 J. Stat. Phys. 83 127 [22] Pathria R K 2008 Statistical Mechanics (Burlington, MA: Elsevier) [23] Leonel E D 2009 Math. Problems Eng. 2009 367921 [24] Leonel E D, Leal da Silva J K and McClintock P V E 2004 Phys. Rev. Lett. 93 014101 [25] Leonel E D, Oliveira D F M and Loskutov A 2009 Chaos 19 033142 [26] Leonel E D and Livorati A L P 2008 Physica A 387 1155 [27] Livorati A L P, Ladeira D G and Leonel E D 2008 Phys. Rev. E 78 056205 [28] Oliveira D F M and Leonel E D 2010 Physica A 389 1009 [29] Loskutov A and Ryabov A B 2002 J. Stat. Phys. 108 995 [30] Lichtenberg A J and Lieberman M A 1992 Regular and Chaotic Dynamics (Appl. Math. Sci. vol 38) (New York: Springer) [31] Bunimovich L A 1991 Chaos 1 187 [32] Eckmann J P and Ruelle D 1985 Rev. Mod. Phys. 57 617 11


J. Phys. A: Math. Theor. 44 (2011) 175102 [33] [34] [35] [36] [37] Manchein C and Beims M W 2009 Chaos Solitons Fractals 39 2041 Manchein C, Rosa J and Beims M W 2009 Physica D 238 1688 Caiani L et al 1997 Phys. Rev. Lett. 79 4361 Firpo M C 1998 Phys. Rev. E 57 6599 Chirikov B V 1979 Phys. Rep. 52 263

A L P Livorati et al

12