Periodic orbits and stability islands in chaotic seas created by separatrix crossings in slow-fast systems

We consider a 2 d.o.f. natural Hamiltonian system with one degree of freedom 
corresponding to fast motion and the other one corresponding to slow motion. The Hamiltonian 
function is the sum of potential and kinetic energies, the kinetic energy being a weighted sum 
of squared momenta. The ratio of time derivatives of slow and fast variables is of order 
$\epsilon $«$ 1$. At frozen values of the slow variables there is a separatrix on the phase 
plane of the fast variables and there is a region in the phase space (the domain of 
separatrix crossings) where the projections of phase points onto the plane of the fast 
variables repeatedly cross the separatrix in the process of evolution of the slow 
variables. Under a certain symmetry condition we prove the existence of many, of order 
$1/\epsilon$, stable periodic trajectories in the domain of the separatrix crossings. Each of 
these trajectories is surrounded by a stability island whose measure is estimated from 
below by a value of order $\epsilon$. Thus, the total measure of the stability islands is 
estimated from below by a value independent of $\epsilon$. We find the location of stable 
periodic trajectories and an asymptotic formula for the number of these trajectories. As 
an example, we consider the problem of motion of a charged particle in the parabolic 
model of magnetic field in the Earth magnetotail.


1.
Introduction. Under rather general assumptions, in the phase space of a mathematical pendulum of a slowly periodically varying length there exist many, of order 1/ε, stable periodic trajectories. Along each of them, the pendulum transits from the oscillatory mode to the rotational one and back again [1,2]. Here ε > 0 is a small parameter characterizing the rate of length variation. In the Poincaré section each of these periodic trajectories is represented by a stationary point surrounded by a stability island whose area is estimated from below by a value of order ε. Hence, the total area of the stability islands is estimated from below by a value independent of ε. Numerical experiments demonstrate that these islands are immersed into a chaotic sea and its measure can be estimated from below by a value independent of ε. However, rigorous estimates of this measure are not known.
The transition from one mode of the pendulum's motion to another one corresponds to a crossing of the separatrix on the pendulum's phase portrait considered at a frozen value of its length. Similar results about stable periodic motions with separatrix crossings are valid in the problem of a particle's motion in a slowly periodically varying double-well potential [1,2].
Analogous results can be obtained for 2 d.o.f. Hamiltonian systems with fast and slow variables (with one degree of freedom corresponding to the fast variables, and another one corresponding to the slow variables). At frozen values of the slow variables we obtain a 1 d.o.f. Hamiltonian system for the fast variables, called the fast system. This system involves values of the slow variables as parameters. Assume that on the phase portrait of the fast system there exists a saddle point, and the separatrices passing through this point form an "8"-like figure. As values of the slow variables change, the projections of the phase points onto the plane of the fast variables may cross the separatrix. Assume that the following symmetry condition is valid: the areas bounded by the two separatrix loops are equal at all possible values of the slow variables. Then, under certain additional conditions, on every energy level there are many, of order 1/ε, stable periodic trajectories crossing the separatrix [3]. Here ε > 0 is a small parameter characterizing the ratio of time derivatives of slow and fast variables. Each one of these periodic trajectories is surrounded by a stability island whose measure is estimated from below by a value of order ε. Hence, the total area of the stability islands is estimated from below by a value independent of ε [3].
The proofs in [1,2] and [3] are based on the study of asymptotic formulas for the corresponding Poincaré maps. In [1,2], these formulas were constructed with the use of asymptotic expressions for the jump of the adiabatic invariant at a separatrix crossing in systems with one and a half d.o.f. [4,5,6] and for the variation of the "angle" variable between separatrix crossings in such systems [7]. In [3], analogous formulas for systems with two d.o.f. [8,9] were used.
In the present paper we consider a natural Hamiltonian system with two degrees of freedom corresponding to fast and slow motions; the Hamiltonian function is a sum of potential and kinetic energies, the kinetic energy is a weighted sum of squared momenta. This system possesses an additional symmetry allowing to locate the periodic trajectories and to find explicitly an asymptotic formula for their number. As main example we consider a model problem [10] of the motion of charged particles in the Earth magnetotail. The results can be also applied to some systems that are not natural but still possess the necessary additional symmetry. In particular, the results are applicable to the model system introduced in [11] to describe the motion of an asteroid near a 3:1 resonance with the Jupiter. The stability islands in this problem were found numerically in [12].
The paper is organized as follows. In Section 2 the equations of motion are presented. These equations contain fast and slow variables. We describe the dynamics of the fast variables at frozen values of the slow variables and the dynamics of the slow variables, averaged over the fast variables. The dynamics is considered on a fixed energy level. Section 3 contains the formulation of the principal analytical results of the paper, Theorems 1 and 2. The assertion of Theorem 1 is that there are many, of order 1/ε, stable in the linear approximation periodic trajectories crossing the separatrix. Explicit asymptotic formulae for the number of these trajectories are given in Section 7. The assertion of Theorem 2 is that the majority of the periodic trajectories whose existence is proved in Theorem 1 are stable and each such stable trajectory is surrounded by invariant tori, which bound a domain of phase volume estimated from below by a value of order ε.
In Section 4 we put forward a standard procedure of approximate description of the dynamics (an improved adiabatic approximation) which describes the behavior of the fast variables on time intervals of length of order 1/ε with accuracy tending to 0 as ε → 0. Such an accuracy is suitable to search for periodic trajectories. In Section 5 we introduce auxiliary variables which are needed to describe the separatrix crossings. These variables are closely related with the variables introduced in Section 4. One of these variables is the improved adiabatic invariant. It is a function of the phase variables which is constant up to terms O(ε 2 ) in the motion far from the separatrix. In Section 6, we present the asymptotic formulae for the Poincaré return map for the problem under consideration. In the adiabatic approximation the behavior of the slow variables is periodic with period proportional to 1/ε. We consider a surface of section to which the phase points return over the time close to this period. This two-dimensional surface of section is parametrized by two variables: the value of improved adiabatic invariant and an angular variable analogous to the phase of the fast motion for phase points near the separatrix. The possibility to construct the return map for such a long period is based on two facts: a) far from separatrix the improved adiabatic invariant changes are O(ε 2 ), while the passage through a narrow neighborhood of separatrices gives rise to a change O(ε); there is a known asymptotic formula for this change, b) the change of the phase of the fast motion with high enough accuracy can be calculated via formulae of improved adiabatic approximation.
In Section 7, we study the problem of existence of stable in linear approximation fixed points of the return map. In Subsection 7.1, we present equations for stable, in linear approximation, fixed points. In Subsection 7.2, it is proved that there are many, ∼ 1/ε , such fixed points, and moreover, asymptotic formulae for the number of such fixed points are presented. The analysis here is based on the reduction of the problem under consideration to the following ergodic problem. On a two-dimensional torus there is a slowly (with velocity O(ε)) moving curve. There is also a point moving on the same torus with slowly varying velocity (the acceleration is O(ε)). One should count the number of times when the point (transversely) crosses the curve on a time interval of length const/ε. A stable, in linear approximation, fixed point of the return map corresponds to each such crossing. Asymptotic formula for the number of crossings can be obtained easily under some generic condition. This completes the proof of Theorem 1. In Subsection 7.3, some symmetry properties of the periodic trajectories are discussed. In Section 8, Theorem 2 is proved.
Section 9 contains an analytical and numerical investigation of an example presenting a model problem of the motion of charged particles in the Earth magnetotail. In Subsections 9.1, 9.2, we describe this model and apply to it the theoretical results of the previous sections. In Subsection 9.3, we present results of the numerical investigation of the problem and demonstrate a good agreement of previous theoretical results with this numerics.
Here (p, q) and (y, ε −1 x) are the pairs of canonically conjugated variables, ε > 0 is a small parameter. All the functions are assumed to be of smoothness C ∞ (but all results are valid for sufficiently high finite smoothness); the functions g and β are positive. The equations of motion have the form: The variables p, q are called fast, while the variables y, x are called slow. The system with Hamiltonian H at frozen values of (y, x) is called the fast system. It contains x as parameter. This is a natural system with one degree of freedom. We assume that the potential U has a unique local maximum at the point q = q s (x) (the subindex "s" stands for "saddle") and that U is a symmetric function of q with respect to q s (x), having two minima, for all the considered values of x.
The phase portrait of the fast system is shown in Figure 1. On it there exist a non-degenerate saddle point C and separatrices l 1 , l 2 passing through C. These separatrices divide the (p, q)-plane into domains G i = G i (x), i = 1, 2, 3. G 1 and G 2 are symmetric w.r.t. the line q = q s (x) and also w.r.t. p = 0. Introduce E by In the domain G 3 the function E is positive, and in G 1,2 (x) it is negative. At the point C and on the separatrices E = 0. Let S = S(x) denote the area of G 1 (x). In the unperturbed system, one can introduce canonical "action-angle" variables (I, ϕ) separately in each G i . The "action" I is a function of p, q, x given by the formula I = I(E, x), where I(e, x) is the area bounded by a trajectory of the fast system, corresponding to E = e, divided by 2π. The action is discontinuous on the separatrix, and hence it is convenient to introduce a modified actionÎ =Î(E, x) coinciding with I at the points of G 1,2 and equal to I/2 at the points of G 3 . On the separatrices and at the point C put I = S(x)/(2π). The functionÎ is continuous.
For an approximate description of the motion in system (2) one can use the adiabatic approximation (see, e.g. [13]). In this approximationÎ = const along a phase trajectory, and the variation of the slow variables y, x is defined by a Hamiltonian system with Hamiltonian εH 0 (I, y, x), where H 0 is the function H expressed in terms of I, y, x. Thus, is the effective potential given by the function 1 2 β(x)p 2 +U (q, x) expressed in terms of I, x. The system with the Hamiltonian H 0 is called the slow system. LetĤ 0 =Ĥ 0 (Î, y, x) denote the Hamiltonian H 0 expressed in terms ofÎ, y, x.
We shall consider trajectories of the slow system on an energy level Assume that on the plane of the slow variables there exists a domain (an annulus) D filled up by closed trajectories of the slow system with H 0 = h 0 . The values ofÎ corresponding to these trajectories fill up a segment Ξ.
Consider on the plane of variables (y, x) the uncertainty curve [11] Γ = {y, x : We assume that every trajectory of the slow system that belongs to D crosses the uncertainty curve at exactly two points. Assume also that at these points Θ(y, x) ≡ gy ∂S/∂x = 0. Note that the points of intersection of a trajectory of the slow system with the uncertainty curve are symmetric with respect to the axis Ox. The uncertainty curve divides D into two subdomains D 3 and D 1,2 . Consider a point (y, x) ∈ D 3 . The corresponding energy level line H = h 0 on the plane of the fast variables (p, q) belongs to the domain G 3 (x) and gives a trajectory of the fast system. For a point (y, x) ∈ D 1,2 the corresponding energy level consists of two trajectories of the fast system, belonging to the domains G 1 and G 2 respectively. For a point (y, x) ∈ Γ the energy level is the union of the separatrices l 1 (x), l 2 (x) and the point C(x). Remark. The value Θ(y, x) is, up to the multiplier ε, the rate of change of the area surrounded by the separatrix loop on the plane of fast variables as the slow variables are evolving according to the slow system. This value plays an important role in the subsequent analysis. In the adiabatic approximation the "action" I remains constant. However, the area surrounded by the separatrix loop is changing. At the points of the uncertainty curve, where Θ(y, x) > 0 (correspondingly, Θ(y, x) < 0), the phase points of the slow system cross the uncertainty curve from D 3 to D 1,2 (correspondingly, from D 1,2 to D 3 ). Accordingly, the projections of the phase points of the exact system onto the plane of the fast variables leave G 3 and enter G 1 or G 2 at Θ(y, x) > 0, and leave G 1 , G 2 and enter G 3 at Θ(y, x) < 0. There is a useful formula (The separatrix l 1 is parametrized by time t of motion in the fast system along it, −∞ < t < +∞. So, on the separatrix ∂E/∂x is a function of t, and the integral in (4) is taken from t = −∞ to t = +∞.) Formula (4) shows that −εΘ(y * , x * ) gives an approximate value of the change of the function E along a segment of a trajectory of the exact system, such that the projection of this segment onto the plane of the fast variables is close to the separatrix l 1 for x ≈ x * , y ≈ y * .
3. Formulation of results. We look for stable periodic trajectories of system (2), belonging to one of the following three families F1, F2, F(1,2). F1. The projection of the phase trajectory onto the plane of the slow variables is a closed curve. Along this curve the value ofÎ oscillates with an amplitude O(ε) around a quantityÎ 0 ∈ Ξ. The period of the motion along this trajectory differs by O(1) from the period of the motion along the trajectoryÎ =Î 0 = const of the slow system. When the projection of the phase point onto the plane of the slow variables is in the domain D 1,2 at a distance larger than c 0 ε from the uncertainty curve, the projection of the phase point onto the plane of the fast variables lies in G 1 . Here c 0 is a large enough positive constant.
F2. The same as F1 with G 1 replaced by G 2 . F(1,2). The same as F1, but the period of motion is O(1)-close to the double of the period of motion along the corresponding trajectory of the slow system, and when the projection of the phase point onto the plane of the slow variables is in the domain D 1,2 at a distance larger than c 0 ε from the uncertainty curve, the projection of the phase point onto the plane of the fast variables is in turn in the domains G 1 and G 2 .
Theorem 1. Under additional generic conditions, on the energy level H = h 0 each of the families F1, F2, and F(1,2) consists of K ε (1 + o(1)) linearly stable periodic trajectories. Among the trajectories of both F1 and F2 there are 1 ε K sm (1+o(1)) trajectories passing through points where y = 0, p = 0 (hence, they are symmetric). Among the trajectories of F(1,2) there are 1 ε K sm (1 + o(1)) trajectories passing through points where y = 0, q = q s (x) + o(1) (hence, they are nearly symmetric; if q s ≡ const, these trajectories pass through the points with y = 0, q = q s and are symmetric). Here K, K sm are positive constants; formulas for them are given in Section 7. The required generic conditions are given in Subsection 7.2 (conditions H1).
This theorem differs from the corresponding result of [3] in the following aspects: here we provide an explicit asymptotic formula for the number of periodic trajectories, we include trajectories of F(1,2) into our consideration, and we describe the locations of the periodic trajectories (defined by the location of a point on the fast plane at y = 0).
The phase space of system (2) is four-dimensional, and the energy level H = h 0 is three-dimensional. Therefore, if there is a two-dimensional invariant torus of the system, it divides the energy level into two invariant domains.
Fix an arbitrary positive constant c 00 . Theorem 2.Under additional generic conditions, among the trajectories of each family, specified in Theorem 1, there are not less than K−c −1 00 ε (and, for symmetric and nearly symmetric trajectories, not less than Ksm−c −1 00 ε ) stable trajectories. Each of these trajectories is surrounded by an invariant torus, which bounds a domain of (3-dimensional) phase volume larger than C −1 1 ε. The variation ofÎ in this domain is smaller than C 2 ε. The required generic conditions are given in Subsection 8 (conditions H2).
Corollary. The total measure of these domains is larger than C −1 3 . Here C i are positive constants, i.e. values independent of ε. This theorem differs from the corresponding result of [3] by the explicit estimate of the number of the stable trajectories.
For any periodic orbit of Theorem 2 there is a Cantor set of tori surrounding it, whose existence is provided by KAM theory. Each one of these tori bounds a solid torus which contains the periodic orbit. The closure of the union of all these solid tori is known as the stability domain associated to the periodic orbit. The intersection of this domain with a surface of a Poincaré section is known as a stability island. 4. Adiabatic and improved adiabatic approximations. In the fast system, the action-angle variables I, ϕ mod 2π are introduced separately in each G i by a canonical transformation of variables. The corresponding generating function W (I, q, x) contains x as a parameter [for brevity, we omit the i subscripts]. We take this function in the form where P is the value of the p-variable along the trajectory with the prescribed value of the action I. In G 1 , G 2 , q 0 (I, x) is the value of q at one of the two points where this trajectory crosses the axis p = 0; in G 1 we take the right one of these points, and in G 2 we take the left one. In G 3 we take q 0 (I, x) = q s (x) and assume that at this point P > 0. In the new variables the Hamiltonian has the form H = H 0 (I, y, x). Now make a canonical transformation of variables (p, q, y, x) → (Ī,φ,ȳ,x) with generating functionȳε −1 x + W (Ī, q, x). The canonically conjugated pairs of variables are (Ī,φ) and (ȳ, ε −1x ). The formulas for the transformation of variables are:φ = ∂W/∂Ī, p = ∂W/∂q,x = x, y =ȳ + ε∂W/∂x.
In the new variables, the Hamiltonian H has the form where In the adiabatic approximation, the dynamics is described by the Hamiltonian H 0 . In this approximationĪ = const along a phase trajectory. One can also construct a canonical, close to the identity, transformation of variables (Ī,φ,ȳ,x) → (J, ψ,ŷ,x) in order to make the terms of order ε in the Hamiltonian independent of the phase (see, for example, [9]). In the new variables, the Hamiltonian takes the form: The term of order ε is absent because it is equal to the average of H 1 with respect toφ, and this average is zero due to the symmetry with respect to the axis p = 0. In the improved adiabatic approximation, the dynamics is described by the Hamiltonian H 0 (J, Y, X). In this approximation J is an integral of motion. With an accuracy O(ε 2 ), the following formula for J is valid (see [8]): The integral here is calculated along a phase trajectory of the fast system passing through the point (p, q); t is the time of motion along this trajectory starting from this point, T is the period of the motion. The function J is the improved adiabatic invariant. In the complete system far from separatrices its value along a phase trajectory is constant with an accuracy O(ε 2 ) on time intervals of order ε −1 . In what follows, it is convenient to use the functionĴ(p, q, y, x) equal to J in the domains G 1 , G 2 and equal to J/2 in G 3 .

5.
Description of separatrix crossing. On the phase plane P sl of the slow variables (y, x), the separatrix is represented by a curve Γ(h 0 ) (see Figure 2). The annulus D filled up with trajectories of the slow system withÎ ∈ Ξ is divided by the uncertainty curve into two domains D 3 and D 1,2 . The ring D intersects the axis y = 0 along two segments S 3 and S 1,2 , which belong to D 3 and D 1,2 respectively. For (y, x) ∈ S 3 the points of the plane of the fast variables on the energy level H = h 0 belong to G 3 , and for (y, x) ∈ S 1,2 the points of the plane of the fast variables on this energy level belong to G 1 or G 2 . Fix an interval Ξ 0 ∈ Ξ, such that the endpoints of Ξ 0 and Ξ are different. Consider the motion in the exact system (2) with initial condition Consider also the motion in the slow system with HamiltonianĤ 0 (Î (0) , y, x) and initial (at This trajectory crosses Γ(h 0 ) at points (x * , y * ) and (x * , −y * ) at slow time moments τ * and τ * * (slow time is τ = εt). We will call τ * and τ * * the slow time moments of the separatrix crossing in the adiabatic approximation.
Assume, for definiteness, that y * > 0, Θ * ≡ Θ(y * , x * ) > 0. Thus, a phase point of the slow system passes from D 3 to D 1,2 at the point (x * , y * ), and at the point (x * , −y * ) it passes from D 1,2 to D 3 . In the exact system, in the process of evolution the projection of the phase point onto the plane of the slow variables (ŷ,x) approaches the curve Γ(h 0 ), and accordingly, its projection onto the plane (p, q) approaches the separatrix. We are interested in the dynamics of phase points that are being captured into G ν , ν = 1, 2 after the separatrix crossing. For such a phase point, the first crossing of Cq-axis in G ν (see Figure 1) occurs near C at time τ = τ * + O(εlnε). Assuming that at the point of this crossing After the separatrix crossing the projection of the phase point onto the plane of the slow variables (ŷ,x) moves towards S 1,2 . When it crosses S 1,2 , the projection of the phase point onto the plane of fast variables is deep inside the region G ν . Denote the value ofĴ at this time moment asĴ (1) .
Then the phase point starts again approaching the separatrix. At τ = τ * * + O(εlnε) the projection of the phase point onto the plane of fast variables crosses Cq-axis in G ν near the point C for the last time before entering G 3 . Assuming that at the point of this crossing After crossing Γ(h 0 ), the projection of the phase point onto the plane (ŷ,x) crosses again the segment S 3 . LetĴ (2) , ψ (2) denote the values ofĴ, ψ at this time moment. Then the projection of the phase point onto the plane of fast variables approaches the separatrix again, crosses it and gets captured into the domain G l , l = 1 or l = 2. Let E be equal to h (2) at the first crossing of the Cq-axis in G l near C. This crossing occurs at time τ = τ * + T 0 + O(εlnε), where T 0 is the slow time period of motion along the trajectoryÎ =Î (0) in the adiabatic approximation. We introduce η (2) = 1 − |h (2) /εΘ * | . 6. The return map. In the energy level H = h 0 the segment S 3 is represented by a piece of a two-dimensional surface {p, q, y, x : H(p, q, y, x) = h 0 , (ŷ,x) ∈ S 3 }. This piece can be parametrized by variablesĴ, ψ. The corresponding Poincaré return map M : (Ĵ (0) , ψ (0) ) → (Ĵ (2) , ψ (2) ) produced by trajectories that pass through G ν is symplectic. Its stable stationary points correspond to stable periodic orbits of the original problem, of period approximately equal to T 0 /ε. It is convenient to study these stationary points using a different set of variables, namely to consider the mapM : (Ĵ (0) , η (0) , ν) → (Ĵ (2) , η (2) , l).
The mapM is the composition of two maps: The results of [8,9] give the following formulas for M (k) . Suppose that Then, as we shall prove immediately, we obtain the mapŝ Here { · } denotes the fractional part, a * = a(x * ). The value of a is given by In (14), (16) one has introduced Here ω ± are the slow time moments when the phase point corresponding to this solution arrives at the separatrix. At these moments S(Y, X) = 2πĴ. (One can take any of such solutions, they differ by a time shift which does not change the value of the integrals in (19), (20)). Conditions (12) , and η (1) , found according to (13), (14) lies in the segment (14) is valid, and analogously for η (2) and expressions (15), (16). All estimates of the form O(·) are uniform on the domain of phase variables under consideration.
Equations (13), (15) follow directly from the formula for the jump of the adiabatic invariant at a separatrix [8] if the assumed symmetry properties hold. Equation (14) follows directly from the formula for the phase change between separatrix crossings [9].
Equations (16)-(18) follow from results of [9]. Indeed, let h + be the value of E when the phase point crosses the axis Cp near C for the first time after exit from G ν . At this time moment the phase point is in G 3 on the positive or the negative part of the axis Cp. Let h − be the value of E when the phase point crosses the same part of Cp for the last time before exit from G 3 . Introduce the notations: According to [9], Section 5) and h + are related as Therefore, Then condition 1/2 < ζ < 1 implies equation (17). Now consider the case l = ν. In this case 0 < ζ < 1/2. We get that ζ and η (2) Therefore, Combining (22) and (23) we obtain equation (16). Condition 0 < ζ < 1/2 implies (18).
We omit the proof of this proposition. The proof consists of repetition of the estimates in [8,9] and additional estimates for derivatives. 7. Linearly stable periodic trajectories.
It follows from assumption H1 that at almost all I ∈ Ξ the winding on the torus with frequency vector γ(I) is ergodic and transversal to Λ(I) at almost all points of Λ(I). Hence, it is possible to calculate, in the main approximation, the number of the points of intersection (cf. [14]). This number is (K + o(1))/ε, where Here l is the natural parameter on Λ(I), and n(l, I) is a unit normal vector on Λ(I).
The factor 1/2 in this formula is due to the fact that the area of the surface of the torus equals 2. The value of the inner integral in (44) equals the product of the length of vector γ(I) and the full length of the projection of the curve Λ(I) onto the normal to vector γ(I).
The curve Λ(I) is the union of two curves Λ nsm (I) and Λ sm constructed according to systems (35)-(37) and (39)-(41) respectively. Hence the value K can also be represented as K = K nsm + K sm , where K nsm and K sm are given by (44) with Λ replaced by Λ nsm (I) or Λ sm . The subscripts "sm" and "nsm" are abbreviations for "symmetric" and "non-symmetric". In the following section we show that indeed, system (39)-(41) gives symmetric or nearly symmetric periodic trajectories.
From the results obtained above we find that if γ 2 γ 1 < 0, K nsm = 0. If γ 2 γ 1 > 0, Λ nsm is a doubled segment of the s 1 -axis on the torus: as η covers the interval (38), this segment is covered twice. Hence, to calculate the length of Λ nsm , one should multiply by four the magnitude of the difference between the values of the function γ 1 αln(2 sin πη) at the middle point and at the end point of the interval (38). Thus, this length equals From (44) we obtain: Now we shall obtain a formula for K sm under the assumptions γ 1 < 0, γ 2 > 0, |γ 1 | > γ 2 ; the results are analogous for the other cases. The projection of a point of the curve Λ sm onto the normal to the vector γ multiplied by the length of γ is The derivative of this function is zero only at the right end of the second interval in (43). Therefore, on both intervals (43) the function f is monotonous, increasing on the first interval and decreasing on the second one. Hence, the product χ(I) of the length of γ and the length of the projection of Λ sm onto the normal to the vector γ equals the sum of the magnitudes of the differences between values of f at the endpoints of the intervals (43): For K sm we have The obtained formulas give the number of periodic trajectories of the family F1 in the main approximation. The number of periodic trajectories of F2 or F(1,2) is the same in the main approximation. For trajectories of F2 the argument is the same. For trajectories of F(1,2) equation (30) should be replaced by and equations (36), (40) should be changed accordingly. The rest of the argument is as before.
If we omit the stability condition | trM ′ | < 2 in the definition of the curves Λ, Λ nsm , Λ sm , the obtained formulas for K, K nsm , K sm give an asymptotic estimate of the total number of periodic trajectories and the number of periodic trajectories with η (1) = η or η (1) = 1 − η satisfying the condition η, η (1) , η (2) . Hence the number of unstable periodic trajectories is bounded from below at least by a quantity of order ε −1 .
In the limit as ε → 0 one can define in a natural way the distribution density of the number of periodic trajectories at a fixed value of I as a function of η. From formulas (44), (35), (36) we find that for trajectories with η (1) = η this density equals const(I)| cot πη|.
Periodic trajectories of the family F(1,2) that correspond to the set Λ sm are nearly symmetric. Using (23) and η + η (1) = 1, we obtain for such trajectories From (50) we find Using (52) we obtain ϕ 0 2π = O(ε 1/3 ln −1/3 ε) mod 1/2. Therefore, possible values of ϕ 0 are O(ε 1/3 ln −1/3 ε) and π + O(ε 1/3 ln −1/3 ε). Hence, at t = 0, when y = 0, the phase point on the plane of the fast variables has a value of q close to q s (x (0) ). If q s ≡ const, the reflection w.r.t. the plane q = q s , y = 0, t = 0 is a canonical transformation. In this case, it can be shown as above that at t = 0 we have q = q s , and the periodic trajectory is symmetric with respect to the plane q = q s . and is surrounded with a Cantor family of invariant curves (see, e.g. [13]). KAMtheory allows also to show that under certain conditions the majority of such fixed points are stable simultaneously, and to give an estimate from below of a size of a stability island around a fixed point which is uniform for the set of these fixed points. To show this we will proceed as follows.
In the new variables the mapM is again a composition of two symplectic maps. It has a fixed pointξ = 0, η = η r . Because of this Denote the principal terms in the r.h.s. of (59) and (60) as g 1 (η r , η r , I r ) and g 2 (η r , η Consider a truncated map: discard terms O(µ) in (61) -(64). Replace in this new map symbols I r , η r withĪ,η and symbol η r replaced with 1 −η. While the original map was defined for a discrete set of parameters I r , η r , the new map is well defined for allĪ ∈ Ξ 0 ,η ∈ [c −1 1 , 1 − c −1 1 ] and has a fixed pointξ = 0, η =η.
Consider the rectangle K = {Ī,η : Let V a be the part of K such that for (Ī,η) ∈ V a the eigenvalues of the fixed pointξ = 0, η =η lie on the unit circle; here a is one of the symbols nsm, sm. On the boundary of V a both eigenvalues equal either +1 or -1. Let λ (a) (Ī,η) be the eigenvalue which has a non-negative imaginary part. The conditions of absence of resonances up to 4th order have the form (λ a ) k = 1, k = 1, 2, 3, 4.
If these conditions are satisfied, then in a neighborhood of the fixed point one can transform the mapM a trunc to the Birkhoff normal form up to terms of the 4th order (see, e.g. [15]). In the coordinates ρ, χ mod 2π which are close to symplectic polar coordinates [15] this form is Here ω a = arg λ a , D a is a smooth function in the domain where the conditions (65) are satisfied. Lower dots denote terms of higher order in ρ. Conditions (65) and condition D a = 0 can be written as one condition of the form F a (Ī,η) = 0, where F a is a smooth function in the domain V a .
We assume that the following condition is satisfied. H2. For allĪ ∈ Ξ 0 we have F a (Ī,η) = 0 for allη such that (Ī,η) ∈ V a except, maybe, at a finite number of points whose location depends onĪ piecewise smoothly; a = nsm, sm. Now we will proceed with the case a = sm, the case a = nsm is completely analogous. The set Λ sm (I) introduced in Section 7.2 is parametrized by values η from a set L sm (I) ⊂ R 1 ((41) implies that L sm (I) consists of two intervals). Fix an arbitrary constant c 2 > 0. Denote asL sm (I) the part of L sm (I) where |F sm (I, η)| > c −1 2 . Condition H2 implies that the setL sm (I) consists of several intervals, and the number of these intervals is bounded uniformly for I ∈ Ξ 0 . The length of L sm (I) \L sm (I) can be made arbitrarily small by taking a sufficiently large value of c 2 . Forη ∈L sm (Ī), the mapM sm trunc in the neighborhood of fixed pointξ = 0, η =η can be written in the form where the estimate O(·) is uniform over the set of valuesĪ,η under consideration. Any symplectic map of the formM sm trunc + O(µ) such that: a)ξ = 0, η =η is its fixed point; b) the estimate O(µ) holds uniformly on the set ofĪ,η under consideration and can be differentiated an arbitrary number of times without changing the order of smallness, can be written in the form This map for small enough µ satisfies KAM stability conditions: there are no resonances up to the 4th order, the coefficient in the normal form that ensures the twist property is different from 0. So, the fixed point is stable. Obviously, the required choice of µ can be made the same for all pairs (Ī,η) under consideration. The original mapM in the neighborhood of the fixed point I r , η r such that η (1) r = 1 − η r + O(µ) (see notations in the beginning of this Section), provided η r ∈L sm (I), can be written in the form (70), (71) withĪ = I r ,η = η r . So, this fixed point is stable for small enough ε, and the required choice of ε can be made the same for all pairs (I r , η r ) under consideration. Now introduceρ = ρ/κ, where κ is a small constant to be defined later. Assuming thatρ = O(1) we will write mapM in the form In the annulus 1 ≤ρ ≤ 2 this map satisfies the conditions of the Moser theorem [16] for the case of a small twist, provided κ is taken small enough and after that ε is taken small enough. So, in this annulus there exists an invariant curve of the map. Obviously, the choice of κ, ε to meet assumptions of Moser theorem can be made the same for all pairs (I r , η r ) under consideration. So, there is an invariant curve at a distance larger than c −1 3 from the pointξ = 0, η = η r . We call the domain surrounded by this curve a stability island of the Poincaré map. The area of this island is larger than c −1 4 . Returning fromξ, η to the original variables I, η, we find that the area of the stability island is estimated from below as c −1 4 ε. The boundary of the stability island (an invariant curve) is a section of the invariant torus on the energy level of the original system. The volume of the domain surrounded by this torus is larger than C −1 1 ε. The variation ofÎ inside this domain is smaller than C 2 ε. This completes the proof of Theorem 2.
9. Example: stable periodic trajectories and stability islands in the problem of the motion of charged particles in the parabolic model of magnetic field.
9.1. The model. The parabolic model (see, e.g. [10]) is one of the simplest models of the magnetic field in the Earth magnetotail. Introduce a right Cartesian coordinate system Oxξq with the origin located in the Earth's equatorial plane, the Earth being located between the origin and the Sun; the Ox-axis lies in the equatorial plane and is directed towards the Earth, and the Oq-axis is orthogonal to the equatorial plane. In this coordinate system the magnetic field vector has the form B = (B 0 q/L, 0, B n ). Here B 0 , B n , L are positive constants. The field lines of this field are parabolas B n x − B 0 q 2 /(2L) = const (see Figure 3).  The vector potential of the magnetic field has the form A = (0, B n x − B 0 q 2 /(2L), 0). The motion of a charged particle in such a field is described by the Hamiltonian system with the following Hamiltonian function: Here m and e are the mass and the charge of the particle, c is the speed of light. The pairs of canonically conjugated variables are (y, x), (η, ξ), (p, q). The Hamiltonian does not depend on ξ, and hence η ≡ const along the motion. Without restricting the generality we put η ≡ 0.
Introduce new pairs of canonically conjugated coordinates (y 1 , ε −1 x 1 ), (p 1 , q 1 ), new time t 1 and a new Hamiltonian H 1 according to the formulas: Here v is a typical velocity of the particle, ρ = cmv/(eB 0 ) is a typical Larmor radius, ε = Bn B0 L/ρ; v, ρ = const. In the new variables the Hamiltonian takes the form (we omit the subscript "1"): Following [10], we consider the case ε ≪ 1 (it is valid for ions). In this case, p, q are the fast variables, and y, x are the slow ones. The phase portrait of the fast system is represented in Figures 4a and 4b for x > 0 and x < 0 respectively. Consider the motion on the energy level H = 1/2. On the configuration plane (x, q) the motion is possible between two parabolas (see Figure 3). The two field lines are also the zero velocity curve (zvc). On the plane of the slow variables the motion is possible in the domain obtained by the union of the half-disk {x, y : x 2 + y 2 ≤ 1, x ≤ 0} and the half-strip {x, y : |y| ≤ 1, x > 0}. The uncertainty curve Γ is the semicircle {x, y : x 2 + y 2 = 1, x ≥ 0}. The phase trajectories on the fast plane that lie inside the domains G 1 and G 2 correspond to the points of the plane of the slow variables with x 2 + y 2 > 1 (see Figure 4a). The phase trajectories that lie in the domain G 3 correspond to the points on the plane of the slow variables with x 2 + y 2 < 1, x > 0.
The phase portrait of the slow system at H = 1/2 is presented in Figure 5 (see [10]). The boundary of the domain of possible motion on the portrait corresponds to the motion with zero value of the "action". For points of the plane of the slow variables with y = ±1, x > 0, the corresponding phase point on the fast plane in the adiabatic approximation always coincides with one of the two stable fixed points of the fast system. For points of the semicircle {x, y : x 2 + y 2 = 1, x < 0}, the corresponding phase point on the fast plane is always at the only stable fixed point of the fast system. The stable fixed point on the phase portrait of the slow system has the coordinates y = 0, x = x e ≈ 0.65. In [10], the chaotic dynamics in the domain filled with trajectories of the slow system that cross the uncertainty curve was found and studied. 9.2. Formulas for the dynamics in the slow system and separatrix crossings. The following formulas are valid for the values S, Θ * , a * , α introduced in Sections 2, 5, and 6 [10]: The modified "action"Î in the domain G i is [10] where K( · ), E( · ) are the complete elliptic integrals of first and second kinds, respectively. It turns out that in the considered specific problem we have: γ 1 < 0, γ 2 > 0, |γ 1 | ≫ 1, γ 2 ≪ 1 (for the definition of γ i , see (26)). Therefore, according to Section 7.2, there are no non-symmetric stable periodic trajectories. Symmetric stable periodic trajectories correspond to the two intervals (43) of the variable η. The corresponding intervals of the variable u = απ cot πη are of equal small length 2/|γ 1 | (see (42)) and in the variable η they are ≈ 2/(απ 2 |γ 1 |) and ≈ αγ 2 2 /(2|γ 1 |), respectively. Hence, the number of stable periodic trajectories is much smaller than the number of unstable periodic trajectories. The minimum of f in (46) is attained at the right endpoint of the second of these two intervals. Therefore, the variation of f on this interval is much smaller than its variation on the first one. Hence, the second interval corresponds to a much smaller number of stable periodic trajectories than the first one. Taking into account the relations between γ 1 and γ 2 , we find for K sm (see (45)) in the main approximation: The quantitiesÎ and x * are related as 2πÎ = S(x * ). Using (76), (77), we obtain Hence, Let v = y 2 * . Then where v 1 and v 2 are the values of v = y 2 * at the endpoints of Ξ 0 . In particular, if we take for Ξ 0 the complete interval of values ofÎ corresponding to the separatrix crossings, we have v 1 = 0, v 2 = 1. Accordingly, for the value K c of the function K sm we find: Here, as usual, B( · ) and Γ( · ) denote the beta and gamma functions respectively. Using the reflection formula for the gamma function [17], we find In the next section we consider solutions with initial conditions p = y = 0, q = q 0 . On the energy level H = 1/2, for the corresponding initial value x 0 of x we have the relation 1/2(x 0 − q 2 0 /2) 2 = 1/2 and, hence x 0 = q 2 0 /2 ± 1. Take the "minus" sign in the latter expression. Calculate for such a solution the value v = y 2 * in the adiabatic approximation. We have 9.3. Results of the numeric studies. For the exact system the 2D surface y = 0 will be used as a surface of the Poincaré section, despite it is not globally transversal to the phase flow. Figure 6 displays a representation of a part of this surface. The horizontal variables are q and p and x is the vertical one. We see that, beyond the fact that it is not a global Poincaré section, to every point (q, p) on the strip |p| < 1 correspond two possible values of x. These effects are illustrated in Figure 7 for ε = 0.04. The same value of ε has been used up to Figure 13. For smaller ε many pictures look quite wild. The initial point ε −1 x = 10, q = 0, y = 0, p > 0 seems to be on an invariant torus. The section by the surface y = 0, disregarding the value of x and the sign ofẏ is shown to the left. The second plot contains only the points for whichẏ > 0. Third one contains points withẏ either positive or negative, but only for the values of x < x e (see Section 9.1). The rightmost plot displays only points withẏ > 0 and x < x e . Using a 3D representation, Figure 8 shows a Poincaré section, disregarding the value of x and the sign ofẏ at y = 0. It has been obtained using initial conditions of the form q = y = 0 and ε −1 x = 3, 4, . . . , 15 and computing 10 4 iterates and then points with the same initial conditions for q and y but with ε −1 x < 3 and step 0.05. We note that most of these points are in chaotic regions and their orbits approach large values of ε −1 x. Typically, the computation is stopped when ε −1 x > 10 4 (or when the number of integration steps using a Taylor method at order 28 exceeds 10 6 steps to compute a Poincaré iterate) to obtain this figure. In the figure there are regions where the density of points is smaller. They correspond to a neighborhood of the lines where y = 0 is no longer transversal to the flow.
Close to x = x e , y = 0 the system has a periodic orbit. It oscillates mildly around the plane y = 0 and intersects it four times. This can be clearly seen in Figure 8, with the corresponding narrow island around it, in the zone of regular motion. Other tiny islands appear inside the chaotic region and the study of them is the main goal of this paper. Some of them can be also seen in Figure 8. However, for a better understanding, we shall consider only the (q, p) projections of the Poincaré section, using the constrains x < x e ,ẏ > 0. This can be seen in Figure 9. The outer part of Figure 9 corresponds to trajectories without separatrix crossings. This part  is mainly filled with segments of invariant curves as it should be by the Arnold perpetual adiabatic invariance theorem. The inner part of Figure 9 is a chaotic sea [10]. The two central holes inside the chaotic region correspond to points coming from large values of ε −1 x (either in q > 0 or in q < 0). The objects inside the chaotic sea which look as segments of fat curves are actually islands with a fixed elliptic point at its center.
A part of Figure 9 is shown in Figure 10, left. Four islands are clearly seen in that plot. A magnification of the rightmost one is shown in Figure 10, right. Because the island is narrow and close to a parabola with horizontal axis, instead of (q, p) we display (q + 0.395p 2 , p).
For a clear identification of the different zones in the configuration plane (ε −1 x, q) we show, in Figure 11, the results presented in Figure 9 in these new variables. The continuous line corresponds to the left boundary of the zero velocity curve (zvc)  Figure 11. Results of Figure 9 projected onto the configuration plane (ε −1 x, q).
As an illustration, Figure 12, top left, shows the first (going from left to right) stable periodic orbit of Figure 10 left, and on the right part the third one. On the bottom we display, for the first of these orbits, the projections on (ε −1 x, y) (left) and on (q, p) (right). On the (ε −1 x, y) plane, it can be seen that the orbit is close to the one theoretically predicted (except for some minor oscillations). On the (q, p) plane the transition of the fast system across the separatrix is clearly seen. We have computed simple periodic orbits (i.e., fixed points of the Poincaré map) with initial conditions on the zvc: y = p = 0. They appear on Figures 9 and 10 on the horizontal axis. Due to the symmetry it is enough to start with q > 0. We note that for y = q = p = 0 there is a periodic orbit fully contained in the (x, y) plane.
In a similar way one can compute periodic orbits starting on the q = 0 axis. Again, due to the symmetry, it is enough to use p > 0. To grasp the wild character of the flow Figure 13 shows, on the left, the first image under the Poincaré map of the line q = 0 with values of p ∈ (−1, 1) and step size 10 −5 . On the right we display, in log 10 scale, the maximum value of ε −1 |x| reached (lower curve) and the return time (upper curve) as a function of p. The plot would be similar taking points of the zvc (p = 0 in this section).
To compute p.o. on the zvc we have restricted our interest to the domain q ∈ [1, √ 2]. (Smaller values of q correspond to extremely large return times.) This is enough to illustrate the theory. This range of q is scanned using a control on the distance between consecutive points, in order to follow carefully all the foldings of the image. Steps in q as small as 10 −12 are accepted. The results of the computation are shown in Table 1. As predicted by the theory, both the total number of p.o. and the number of stable ones (s.p.o.) behave as ε −1 .
The p.o. in p = 0 appear also (if we display the period or the maximal value of x along the orbit) in bumps similar to the ones displayed in Figure 13, right. Decreasing ε the number of bumps increases like ε −1 . On each bump the behavior of the system is similar.   A representation of the location of the p.o. is displayed on Figure 14, left. On it the cumulative number of p.o. which are found from q = 1 on, multiplied by the corresponding value of ε, is displayed as a function of q. The curves are shown for ε = 0.01, 0.0025, 0.001 and 0.00025. The smaller is the value of ε, the smaller is the effect of the bumps. Then, on the right part of Figure 14 we have collected the distribution of p.o. inside a bump as follows. If ε is small the bumps are also small. The extrema of the bump can be taken as the places where the return time has a local minima (or, equivalently, the maximal x is minimum or, also, the trace of the nearby p.o. is also a minimum). If q 1 and q 2 are the extrema of some bump (q 1 < q 2 ), and q * denotes the position of a p.o., it can be referred to the bump by means of the parameterq = (q * − q 1 )/(q 2 − q 1 ). Collecting all the values ofq for the different p.o. we can get the desired distribution. This is plotted to the right in Figure 14. Despite a total of 359 complete bumps is found for q ∈ [1, √ 2] for ε = 0.00025 (this reduces to 89 for ε = 0.001), only the bumps for q ∈ [1, 1.1] have been used. The reason is that for that range of q the computed value ofq is close enough to the theoretical parameter η introduced in Section 5. In this range there are 58 full bumps which contain 10373 p.o. We selected 200 intervals of equal length in theq variable. The scale has been normalized to have total area equal to 1. In the same plot, and in discontinuous lines, the theoretical density of the p.o. as a function of η is also shown. It is given by an expression of the form (49) where c = 2(γ 1 + γ 2 )/(πγ 1 γ 2 ) (in notations of Section 7) and C is a normalization constant. The approximate values c = 1.8, C = 0.3 have been used. Due to the relatively small number of p.o. the experimental density shows some irregularities, but the global agreement is remarkable. We recall that close to the endpoints of the interval the theoretical approximation fails. The density becomes zero at 1 π tan −1 (1/c) ≈ 0.16, in good agreement with numerical results.
Concerning the s.p.o., in each bump one should find two narrow intervals in η (43) which can contain such orbits. The largest one is close to η = 1/2. However, as the intervals are narrow and the transformation fromq to η depends also on the location of the bump, it is not suitable to look at results considering all the bumps together. As before, it is enough to use only a part of the bumps (for instance, with q 0 ∈ [1, 1.1]). It is checked that the stable p.o. are located in the interval q ∈ [0.5, 0.515]. A very narrow interval of s.p.o. appears also nearq = 0.2, but it contains less than 2% of the total number of s.p.o. This is in agreement with the theoretical prediction in Section 9.2.
We compare the total number of s.p.o. to the theoretical predictions as follows. Orbits starting on the zvc for q 0 ∈ [1, √ 2] cross the separatrix (x 2 + y 2 = 1) at values of y ∈ [y 1 , y 2 ]. Let v i = y 2 i . According to (84), we have in the adiabatic approximation Calculating the integral in (47) with these values of v 1,2 we obtain K sm ≈ 0.099. The quantity K sm /ε is, in the main approximation, the number of s.p.o. of the family F1 (or F2) that start from the zvc with q 0 ∈ [1, √ 2] and q 0 ∈ [− √ 2, −1]. It is also equal in the main approximation to the number of s.p.o. of the two families that start from the zvc with q 0 ∈ [1, √ 2]. For K sm /ε we find values which are in good agreement with the numerical results (see Table 1).
Appendix. The planar periodic orbit. Here we discuss the special role played by the periodic orbit sitting on (q, p) = (0, 0). As it is contained in the (x, y) plane, we shall call it the planar p.o. If y(0) = 0, it is given by x = cos(εt), and the normal variational equation is given (as a second order equation) byξ = cos(εt)ξ. This is a particular case of Mathieu equation, with zero average. Let M = π 0 sin(τ )dτ ≈ 2.39628. For small ε the monodromy matrix can be approximately written as the product of rotation by the angle M/ε and an hyperbolic matrix of dominant eigenvalue exp(M/ε). Figure 15, left, displays the variation of the trace tr of the monodromy matrix as a function of 1/ε. In the horizontal axis we display 1/ε, while in the vertical one sinh −1 (tr) is shown. It is straightforward to show that, despite infinitely many intervals of stability occur they are extremely narrow. More precisely: stability domains appear near values of 1/ε of the form M −1 (k + 1/2)π, k ∈ Z, k > 0, while the width of these intervals is, up to a constant factor, of the form exp(−M/ε). In each one of these stability intervals there is an invariant island around (0, 0) in the (q, p) plane. As an example, Figure 15 right shows the island for ε = 0.5155, in the stability domain with the left endpoint at ε = 0.5136759883096067 and width ≈ 2.25 × 10 −3 . Even for that large value of ε the width of the island in q is small (less than 0.002). Furthermore there exists a period 4 p.o. (displayed as black big points) which is hyperbolic. The presence of nearby hyperbolic periodic orbits has been checked for other values of ε in the stability domains. Hence, when the planar p.o. becomes stable (for these narrow intervals) there appear other unstable p.o. which play a similar role. Using an index argument in a suitable ball containing the island inside, one can prove that the global dynamics outside this ball is only affected for very narrow intervals of initial conditions. From now on we shall not consider these stability intervals.
Dependence of the return time on the initial value of q contains certain characteristic structures that we call bumps, see Figure 13, right. Figure 16 shows, in the upper part, a sample of results, and several details in the lower part. Concretely, on the left upper plot we present the results for trajectories with initial conditions at p = 0 and q ∈ [0, √ 2] for ε = 0.3, 0.2, 0.1, 0.04. We plot the logarithm of the return time as a function of the initial value of q. The smaller ε is, the higher are the bumps and the smaller is the distance between endpoints of a bump, as it was said before. For instance, for ε = 0.04 the maximal return time from the zvc exceeds 10 7 . In the central upper plot we show, for q ∈ [1, √ 2] the logarithms of the return time multiplied by ε for ε = 0.04, 0.01, 0.0025. The curves on these plots tend to a limit curve when ε → 0. On the right upper part we show the end of the big bump close to q = 1.15 for ε = 0.04, displaying the return time as a function of q on the narrow window [1.15027, 1.1502744]. Five ends of small bumps can be seen to the left of the endpoint of the big bump, one of them very close to this endpoint. All this detailed structure is not visible in the previous plot (the middle upper one). In the lower part of the figure, and from left to right we show a detail of the end of the leftmost small bump of the previous plot (window [1.150270992 ± 2 × 10 −9 ]), a view of the end of the rightmost small bump together with the end of the big bump (window [1.150274085 ± 5 × 10 −9 ]) and a great enlargement of the end of the big bump (window [1.150274092745 ± 25 × 10 −12 ]. The vertical jumps in the return time, clearly visible in last plot, are due to the fact that p x = 0 is not everywhere transversal to the flow, as it was discussed (see Figure 7 right, for instance). To the right of the endpoint of the big bump the pattern of small bumps is almost symmetric to the one displayed here. Now we pass to the reasons for the creation of big and small bumps in geometric terms. Taking, for instance, an initial point on the zvc, it can approach the separatrix x 2 + y 2 = 1 moving close to the planar p.o. Then it can move back to smaller values of x, or it can enter one of the branches of the region x > 1, either the upper one or the lower one. A critical orbit would tend to the planar one along its stable manifold. Hence, looking at the maximal value of x along the orbit, it would never exceed 1. Furthermore, moving on the stable manifold of the planar p.o. it would reach a point with y = 0,ẏ > 0 close to x = −1.
Hence, the endpoints of the big bumps can be seen as the place where there is a minimum of the maximum value of x along the orbit (as a function of the initial value of q in the zvc) or as the place where orbits switch from entering the upper branch to the lower one or vice versa (see below for an additional explanation). Or, in other words, to switching between entering one or the other of the lobes at the crossing of the separatrix.
Small bumps are related to the crossing of the separatrix when we go outside, to the circulation region, that is, x decreases and enters the domain x < 1. An orbit entering from one of the branches to the region x < 1 can also approach the planar p.o. along the stable manifold. This occurs depending on the phase it has when it reaches the maximal value of x in one of the branches. As it changes quickly for ε small, there are many small bumps inside one big bump. Again, entering x < 1 we can pass close to the planar p.o. and leave its vicinity following one of the two branches of their unstable manifold.