Geometry of complex instability and escape in four-dimensional symplectic maps
GGeometry of complex instability and escape in four-dimensional symplectic maps
Jonas St¨ober
1, 2 and Arnd B¨acker
1, 2 Technische Universit¨at Dresden, Institut f¨ur Theoretische Physik and Center for Dynamics, 01062 Dresden, Germany Max-Planck-Institut f¨ur Physik komplexer Systeme, N¨othnitzer Straße 38, 01187 Dresden, Germany (Dated: September 3, 2020)In four-dimensional symplectic maps complex instability of periodic orbits is possible, whichcannot occur in the two-dimensional case. We investigate the transition from stable to complexunstable dynamics of a fixed point under parameter variation. The change in the geometry ofregular structures is visualized using phase-space slices and in frequency space using the exampleof two coupled standard maps. The chaotic dynamics is studied using escape time plots and bycomputations of the invariant manifolds associated with the complex unstable fixed point. Basedon a normal-form description, we investigate the underlying transport mechanism by visualizing theescape paths and the long-time confinement in the surrounding of the complex unstable fixed point.We find that the escape is governed by the transport along the unstable manifold across invariantplanes of the normal-form. I. INTRODUCTION
There are different ways in which orbits of a dynamicalsystem may become unstable under variation of some pa-rameter. One famous example is the Hamiltonian-Hopfbifurcation as has first been studied for the triangularequilibrium points of the planar circular restricted three-body problem [1, 2], for which instability occurs beyonda critical mass ratio [3]. This is also found for manyother examples in celestial and galactic dynamics [4–12]for the hydrogen atom [13–15], in the context of molecu-lar dynamics [16, 17], and is also of relevance to particleaccelerators [18]. The impact of the Hamiltonian-Hopfbifurcation on the phase space geometry has been stud-ied in much detail in Refs. [19–21]. Additional insightis provided by computations of invariant manifolds andnormal-form descriptions [22–25]. For further results seee.g. Refs. [26–29]. The impact in quantum mechanicalsystems has been investigated in Ref. [30].Often its is helpful to reduce the time-continuousdynamics to a discrete-time mapping by means of aPoincar´e section. For conservative Hamiltonian systemswith three degrees of freedom this leads to the study offour-dimensional ( ) symplectic maps, which are there-fore of importance of many areas of physics. Similar tothe Hamiltonian case, a transition from stable to com-plex unstable dynamics is possible for (and higher-dimensional) symplectic maps [31, 32]. This has beeninvestigated in detail in the pioneering work [33, 34] fora variant of the coupled standard map [35]. In such atransition to complex unstable dynamics two eigenvaluepairs of the linearized dynamics collide on the unit circleand afterwards form a so-called Krein quartet. This mayonly happen, if the Krein signature is mixed [31]. A dis-tinctive feature is the spiraling motion in the surround-ing of a complex unstable periodic point [6, 30]. More-over, it was found that commonly an extended regionaround a complex unstable fixed point emerges to whichthe dynamics is confined for rather long times [11, 36–38]. Important approaches to understand the complex unstable dynamics are based on computations of the in-variant manifolds [36, 38, 39] and normal form descrip-tions [15, 40, 41]. Hamiltonian-Hopf bifurcations havealso been studied in much detail for reversible maps, seee.g. Refs. [42, 43].In this paper, we investigate how the transition fromstability to complex instability of a fixed point affectsthe geometry of invariant objects in its surrounding inthe phase space of a symplectic map. This transitionis accompanied by the possibility that orbits can escapefrom the vicinity of the fixed point which is quantified bythe average escape times of an ensemble of orbits. Theunderlying escape mechanism is investigated in terms ofthe geometry of the stable and unstable manifolds. Weprovide evidence that the escape occurs across the invari-ant planes of the normal-form description showing thatit is a genuinely higher-dimensional mechanism.The text is organized as follows. In Sec. II we recallsome fundamental properties of linear stability of fixedpoints and the requirements for complex instability in symplectic maps. Section II C summarizes a normal-formdescription for the transition to complex instability as in-troduced in Ref. [40]. In Sec. III we introduce a variantof the four-dimensional coupled standard map and definea set of parameters for investigating the transition fromelliptic-elliptic stability to complex instability. We visu-alize the dynamics in the phase space using phase-space slices [44] which is complemented by a frequencyspace representation [45–47]. The escape dynamics is in-vestigated in Sec. IV for an ensemble of initial conditionsclose to the complex unstable fixed point. To explainthe underlying mechanism we compute the stable andunstable manifolds associated with the complex unsta-ble fixed point by utilizing the parametrization method[48–50]. The dynamics of the ensemble suggests that theescape occurs across invariant planes of the correspond-ing normal-form description. Section V gives a summaryand outlook. a r X i v : . [ n li n . C D ] S e p II. COMPLEX UNSTABLE DYNAMICSA. Linear stability in 4D maps
In this section we collect some important results on thestability of fixed points in symplectic maps [31], theKrein collision [31, 51–53] and its normal-form descrip-tion [40, 54, 55]. A map M : R → R is called symplec-tic if its Jacobian matrix D M fulfills D M T J D M = J ,where J = (cid:18) − II (cid:19) is the 4 × I being the 2 × M ) = 1. The dynamics in the vicinity of a fixedpoint, i.e., a point z ∗ that satisfies M z ∗ = z ∗ , is given bythe linearized map D M . The symplecticity of M impliesthat the characteristic polynomial P ( λ ) of D M is reflex-ive so that coefficients of P come in palindromic form.For a symplectic map this can be written as P ( λ ) = λ − Aλ + Bλ − Aλ + 1 , (1)where A = tr (D M ) and 2 B = A − tr (D M ). As conse-quence, the eigenvalues λ j with j ∈ { , } are restrictedto either hyperbolic pairs λ j , λ − j ∈ R , elliptic pairs of λ j , ¯ λ j ∈ C with | λ j | = 1 or a Krein quadruplet of com-plex eigenvalues λ, λ − , ¯ λ, ¯ λ − ∈ C with | λ | (cid:54) = 1.This gives a total of four possible stability types,namely elliptic-elliptic (EE), elliptic-hyperbolic (EH),hyperbolic-hyperbolic (HH) and complex instability(CU). These stability types can be distinguished by in-troducing the stability index of an eigenvalue pair ρ = λ j + λ − j and reducing the characteristic polynomial inEq. (1) to R ( ρ ) = P ( λ ) λ − = ρ − Aρ + B − . (2)As shown in Ref. [31], different regimes of stability followfrom Eq. (2) in dependence on A and B . The linearizedmap D M is spectrally stable if and only if all roots of R ( ρ ) are real and within the interval [ − , R ( ±
2) = 0 yields two stability boundaries, namely B = ± A − . (3)Crossing either of these boundaries corresponds to asaddle-center (SC) or a period-doubling (PD) bifurca-tion, respectively. Another boundary corresponds to theroots of R ( ρ ) becoming complex, which occurs whenthe discriminant of the reduced characteristic polynomial∆ R ( ρ ) = ( ρ − ρ ) = 0. This gives the so-called Kreinparabola (KP) B = A / . (4)The possible stability types for an arbitrary fixed pointof a map in dependence of A and B can be dis-played in the so-called Broucke diagram [1, 31], see Fig. 1.The three stability boundaries SC, PD, and KP lead to − − − − II HH
IH EHIE EECU S C P D K P AB FIG. 1. Stability of a fixed point in dependence on the coef-ficients A and B of the characteristic polynomial (1) of thelinearized map D M . The regions correspond to combina-tions of elliptic (E), hyperbolic (H), and inverse hyperbolic(I), or complex unstable (CU). The regions are seperated bythe period-doubling line (PD), saddle-center line (SC), andthe Krein parabola (KP). seven stability regions corresponding to complex insta-bility (CU) and the different combinations of the elliptic(E), the hyperbolic (H) case, and the inverse hyperbolic(I) case, for which the eigenvalue pair lies on the negativereal axis. The corresponding arrangement of the eigen-values of the linearized map are shown as small insets.For an EE fixed point the surrounding consists of atwo-parameter (Cantor) family of tori as expectedfrom Kolmogorov-Arnold-Moser (KAM) theory. The tori are organized around one-parameter (Cantor) fami-lies of elliptic tori. These families are commonly re-ferred to as Lyapunov families, based on the analogy tothe Lyapunov center theorem for Hamiltonian flows [56].Such families of tori have been studied in detail, seee.g. Refs. [36, 57–62]. As the families of elliptic toriform the ‘skeleton’ of the surrounding regular dynamics,they allow for a convenient way to understand the changein geometry occurring when an EE fixed point becomesCU, as will be illustrated below in Sec. III B. B. Krein collision
As seen from Broucke’s diagram in Fig. 1, there areonly three possible ways to enter the CU regime, namelythe transition from a) the elliptic-elliptic (EE), or b)the hyperbolic-hyperbolic (HH or II) stability regionsthrough the Krein parabola, or c) through the intersec-tion points of the Krein parabola with either the saddle-center or the period-doubling boundary at (
A, B ) =( ± , α , which con-trols the transition. For α >
0, two elliptic eigenvaluepairs approach each other on the complex unit circle un-til they coalesce at α = 0. For α <
0, the eigenvaluessplit off the unit circle and form a Krein quadruplet.Whether the eigenvalue pairs of an EE fixed point for agiven map can leave the unit circle or pass through eachother while staying on the unit circle depends on the so-called Krein signature. This is given by the signature( m + , m − ) of the quadratic form q ( x ) = x T J D M x, (5)which can for example be computed numericallyfrom the eigenvalues of the symmetric matrix (cid:0) J D M + ( J D M ) T (cid:1) , where m + is the number ofpositive and m − is the number of negative eigenvalues.If m + = 0 or m − = 0 then the fixed point cannotloose stability and stays elliptic-elliptic. Conversely,the fixed point may loose its stability and becomecomplex unstable if the signature is mixed. Note thatthe quadratic form Eq. (5) allows the construction of aninvariant of the linearized dynamics as q ( x ) = x T J D M x = (D M x ) T J D M (D M x ) (6)is preserved under D M [33].The geometric interpretation of the Krein signaturebecomes more clear when considering the signature of amultiplier λ on the unit circle, σ ( λ ) = sgn q ( u ) , (7)where u is any real vector in the eigenspace of λ . Ifeigenvalues with the same signature collide on the unitcircle, they cannot split off to form a Krein quartet. Ex-plicitly, consider a symplectic map which is uncou-pled, i.e. M ( p , p , q , q ) = ( p (cid:48) , p (cid:48) , q (cid:48) , q (cid:48) ) with ( p (cid:48) , q (cid:48) ) = M ( p , q ) and ( p (cid:48) , q (cid:48) ) = M ( p , q ). Then using thequadratic form (5) and (1 , , ,
0) and (0 , , ,
0) as vec-tors of the corresponding eigenspaces the signatures aregiven by σ ( λ i ) = sgn((D M i ) ). Therefore the fixedpoint can only become complex unstable under somegeneric coupling if [63, 64]sgn((D M ) ) sgn((D M ) ) < . (8)This reflects the counter-rotating nature of the dynamicsin the two independent subspaces, similar to the Cherry-Hamiltonian describing two counter-rotating harmonicoscillators [65].Furthermore, a mixed Krein signature implies thatthe linearized map of the coalesced eigenvalues has non-trivial Jordan blocks of the shape m + × m + and m − × m − while the matrix can be diagonalized if the signature α > ϕ Re λ Im λ α = 0(b) ϕ Re λ Im λ α < ϕ Re λ Im λ FIG. 2. Krein collision of two elliptic eigenvalue pairs (red andblue circles) in dependence of α . The eigenvalues coalesce for α = 0 and split off the complex unit circle for α < ϕ thelocation of the Krein collision is move along the unit circle. is positive or negative definite. Thus, the linearizationtakes either the form [33] λ λ λ
10 0 0 ¯ λ or λ λ λ
00 0 0 ¯ λ (9)where λ = e i θ and θ ∈ ]0 , π [. Beside this, in case b) thesignature is always mixed. Thus for an II or HH fixedpoint there is no constraint to enter the CU region. C. Normal form description
To understand the geometry of regular and invariantstructures around a CU fixed point, it is helpful to con-sider a non-linear normal-form description [40], of whichwe now summarize the main aspects. Consider a sym-plectic map M , x (cid:48) = M ( x , α, ϕ ) , (10)with x , x (cid:48) ∈ R and parameters α, ϕ ∈ R . The fixedpoint is assumed to be at the origin x = such that M ( ; α, ϕ ) = 0 for arbitrary α and ϕ . Furthermore,the eigenvalues of the linearized map D M ( ; 0 ,
0) are as-sumed to coalesce ot λ = exp ( ± i θ ) with θ = 2 πν andirrational ν ∈ ]0 , / [. Note that the case of the ratio-nal Krein collision is for example considered in Ref. [55].The collision is controlled by the parameters α and ϕ asshown in Fig. 2. The parameter α controls the transitionfrom the elliptic-elliptic eigenvalue pair for α > α <
0. The angle ϕ ro-tates the angle of the Krein collision on the complex unitcircle.In case of the irrational Krein collision with α = 0 and ϕ = 0, the linearized map has non-trivial Jordan blocksand can be brought into a Williamson normal-form L by a symplectic transformation TT − D M T = L ( ; 0 ,
0) = (cid:18) R θ (cid:15)R θ R θ (cid:19) (11)where (cid:15) = ± R θ = (cid:18) cos θ sin θ − sin θ cos θ (cid:19) . (12)For α (cid:54) = 0 and ϕ (cid:54) = 0, the Williamson normal-form has atransversal two-parameter unfolding, i.e., there is a two-parameter family of matrices that preserve the symplec-tic form and describes the transition from stability tocomplex instability via the Krein collision given by [55] L = L ( , α, ϕ ) = (cid:18) (1 − (cid:15)α ) R θ + ϕ (cid:15)R θ + ϕ − αR θ + ϕ R θ + ϕ (cid:19) . (13)With that, the transformed map (cid:102) M in the new coordi-nates y can be represented as a formal power series y (cid:48) = (cid:102) M ( y , α, ϕ ) ≈ L + Φ ( y , α, ϕ ) + . . . (14)where Φ j ( y , α, ϕ ) are vector-valued polynomials of degree j . In Refs. [40, 54, 55] it is shown that Eq. (14) canbe normalized by utilizing a symplectic diffeomorphismΨ j ( y ) such that Ψ − j ◦ (cid:102) M ◦ Ψ j is in normal-form withrespect to L up to order j for arbitrary j ∈ N .As a result, one gets the non-linear normal-form [40] (cid:18) x (cid:48) y (cid:48) (cid:19) = (cid:20) (1 − (cid:15)h ) R θ + ν (cid:15)R θ + ν − hR θ + ν R θ + ν (cid:21) (cid:18) xy (cid:19) (15)with ( x, y ) = ( x , x , y , y ) ∈ R . The parameters h = α + b X + b I + . . .ν = ϕ + b X + b I + . . . . are derivatives of a deduced Hamiltonian generating func-tion with respect to the coordinates X = x + x and I = y x − x y , respectively. For our purposes we trun-cate the series of the generating function after quadraticorder and obtain (cid:101) ν = ϕ and (cid:101) h ≈ α + bX where (cid:101) h is scaledwith respect to X such that b = ± S -equivariant where the symmetry transfor-mation acts as rotation on all coordinates ( R γ x, R γ y ) for γ ∈ [0 , π [. The corresponding invariant of Eq. (15) is I ( x (cid:48) , y (cid:48) ) = I ( x, y ) = y x − x y . Consequently, the nonlinear normal form map can be reduced further byintroducing new coordinates.Hence, we take advantage of the symmetry and visu-alize the dynamics of Eq. (15) in the hyperplane x = 0,see Fig. 3. Note that the full dynamics can be re-obtained by applying the symmetry operation, i.e. bysimultaneous rotation in the x and y coordinates, seeRef. [40, Eq. (3.1)]. For the sake of clarity, we stickto the half-space with x ≥ x , y ) (cid:55)→ ( − x , y ).Furthermore, without loss of generality we fix the pa-rameters (cid:15) = 1 and b = 1. Firstly, we consider the case x y y I = 0 I < x y y I = 0 I <
FIG. 3. The reduced Poisson map from Eq. (15) in ( x , y , y )coordinates. The sphere in the origin denotes the trivial fixedpoint while the gray and the blue planes visualize the I = − .
015 and the I = 0 plane, respectively. The shown orbitscorrespond to the same plane as their color indicates. Thenon-trivial periodic points of the reduced map are depictedas orange and yellow dots for the EE case (a) for α > α < I = 0, i.e. 0 = I = − x y . Without loss of general-ity, we choose y = 0 and Eq. (15) reduces to a map f ( x , y ) (cid:55)→ ( x (cid:48) , y (cid:48) ) x (cid:48) = | ( gx + y | (16a) y (cid:48) = ( (cid:101) hx − y ) sign( gx + y ) (16b)with g = 1 − (cid:15)h . This map has two periodic points,namely a trivial fixed point at (0 ,
0) which is the originalfixed point of M and for α < (cid:112) − α/b, α as expected. In contrast, the non-trivial peri-odic point only exists when α ≤ I = 0 plane correspondsto the typical behavior of a period-doubling bifurcationin a symplectic map, for which a periodic point loosesits stability and a stable periodic point of the twice theperiod is created, see e.g. [66–69].For the second case I (cid:54) = 0, the coordinate y is given bythe invariant I . Thus, Eq. (15) reduces to a map withall structures living on a hypercolic cylinder y = − I/x in the reduced phase space. The map takes the form x (cid:48) = (cid:113) ( gx + y ) + I / x (17a) y (cid:48) = ( gx + y )( y − (cid:101) hx ) + I / x x (cid:48) . (17b)In this case, there is only one non-trivial period-twopoint, which is given by an implicit equation that wesolve numerically.Figure 3(a) shows the reduced phase space in( x , y , y ) coordinates for α >
0, i.e., the stable case.The red sphere represents the trivial fixed point which iselliptic-elliptic in this case. The blue and the gray planesas well as the orbits in the same color correspond to I = 0and I = − . α is positive,there exits only one fixed point in the I = 0 plane. For I >
I < α < I (cid:54) = 0plane are detached from the origin similar to a period-doubling bifurcation in a map. In this way, this familywith its surrounding stable tori forms a foliated tube-like object in phase space. III. TRANSITION TO COMPLEX INSTABILITYA. 4D map with CU fixed point
The usual standard map [70, 71], which has beeninvestigated in much detail, see e.g. Refs. [37, 44, 72–74], does not allow for CU fixed points. A modified standard map has been introduced in Ref. [33], whichis inspired by the Cherry-Hamiltonian describing twocounter-rotating harmonic oscillators [65]. As exemplarysystem to study the transition from EE to CU stabilitywe use a variant of such two coupled counter-rotating standard maps given by the map M ( p , p , q , q ) (cid:55)→ ( p (cid:48) , p (cid:48) , q (cid:48) , q (cid:48) ) as p (cid:48) = p + K π sin 2 π ( q (cid:48) ) + K π sin 2 π ( q (cid:48) + q (cid:48) ) (18a) p (cid:48) = p + K π sin 2 π ( q (cid:48) ) + K π sin 2 π ( q (cid:48) + q (cid:48) ) (18b) q (cid:48) = q + p (18c) q (cid:48) = q − p , (18d) where K and K are the kicking strengths of the two subsystems and K determines the coupling betweenthem. The phase space is restricted on the torus, i.e.( p , p , q , q ) ∈ [ − . , . × [0 , with periodic bound-ary conditions. Note that the counter-rotating charac-ter of the two uncoupled subsystems in ( p , q ) and( p , q ) is due in the negative sign of the second momen-tum p in Eq. (18d) which ensures that condition (8) isfulfilled. This sign is the only difference to the usual standard map, as introduced in Refs. [70, 71]. This maphas also been investigated in Ref. [75], though with thenegative sign in Eq. (18c) instead of Eq. (18d).We will focus on the central fixed point at z ∗ =(0 , , / , / ) in the following. Its stability coefficientsare A = − K + K + 4 , (19a) B = − K K + K K − K + K K + 2 K + 6 . (19b)Figure 4 shows the stability diagram for fixed coupling K = 0 . K and K . The fixed point iscomplex unstable in the region between the two straightlines K = − K and K = 4 K − K . (20)The saddle-center and the period-doubling boundaries,Eq. (3), lead to the hyperbolae − KK K − K and − KK + 4 K − K − K + 4 . (21)In order to investigate the transition from EE to CUstability, we choose the EE region with positive kickingparameters and fix K = 0 . K is varied. Thesix equidistant parameters K = 0 . , . , ..., . K, K , K ) = (0 . , . , . λ, λ − , ¯ λ, ¯ λ − ) of D M where λ = exp ( β + i θ ) with β ∈ R + and θ ∈ [0 , π [, see Sec. II A. The correspondingeigenvectors ( ξ , ξ , ¯ ξ , ¯ ξ ) can be written as ξ j = u j + i v j with u j , v j ∈ R and j = 1 ,
2. The stable and unstableinvariant subspaces of the linearized map are spanned by u , v and u , v , respectively.From this one key feature of the dynamics in the sur-rounding of a CU fixed point follows: Under the lin-earized dynamics these eigenvectors evolve as ξ nj = λ nj ξ j and consequently provides the evolution in the stable andunstable subspaces by [21] u ( n ) j = exp ( ± βn ) (cos ( nθ ) u j − sin ( nθ ) v j ) (22a) v ( n ) j = exp ( ± βn ) (sin ( nθ ) u j + cos ( nθ ) v j ) , (22b)where the positive sign corresponds to j = 1 and thenegative to j = 2. Any point z in the phase spacecan be expressed in the basis of the eigenvectors, i.e. z = − . . . . − . − . . . . . . K K (A)(F) FIG. 4. Stability of the fixed point (0 , , / , / ) for fixed K = 0 . K and K . The magnificationshows the selected parameters for the transition from EE toCU, (A) K = 0 .
31, (B) K = 0 . K = 0 . c u + c v + c u + c v with coefficients c , c , c , c ∈ R .These coefficients can be determined with the help of thebasis of the dual space of the matrix of eigenvectors [44].Using the time evolution of the eigenvectors Eq. (22) al-lows for obtaining the linearized dynamics of an orbitfor a given initial condition. Apparently, the underlyingdynamics is governed by an expanding/contracting partand a rotating part which leads to a spiraling motion asillustrated in Fig. 5. If c or c are different from zero, theexpanding dynamics will asymptotically dominate. Notethat this provides a good description for some limitednumber of iterations of the map M only, beyond whichthe nonlinear dynamics becomes relevant, as can be seenby the deviations between the real orbit depicted as col-ored spheres and the linerized dynamics shown as blackcurve in Fig. 5. minmax p FIG. 5. Spiraling motion of an orbit started close to theCU fixed point. Shown are the ( p , q , q ) coordinates with p encoded in color of the first 170 iterates of the point(0 , , . , .
5) + µ for µ = 10 − . The initial spiraling motionis well described using the linearized dynamics, Eq. (22), asshown by the black curve. From the 160th iterate deviationsbecome visible in the plot. B. 3D phase-space slice
To get an intuition for the dynamics of the transitionfrom EE to CU stability of the fixed point in phase space,we use a phase-space slice [44]. The idea is to reducethe phase space by one dimension by considering a hyperplane Γ and determining those points of an orbitthat fulfill the slice conditionΓ ε = (cid:8) ( p , p , q , q ) (cid:12)(cid:12) | p | ≤ ε (cid:9) . (23)For the resulting points the coordinates ( p , q , q ) aredisplayed in a plot. The parameter ε , i.e., the thick-ness of the slice, controls the resolution. Smaller valuesof ε require longer orbits to obtain the same number ofpoints in the slice as the slice condition (23) is fulfilledless often. For all phase-space slice plots in this paperwe choose ε = 10 − . Typically f -dimensional objects inthe full phase space appear as ( f − phase-space slice. For example tori leadto two (or more) separate (but dynamically connected)rings in the phase-space slice and tori lead to two(or more) points in the slice. For further examples, alsoincluding more general slice conditions, and detailed dis-cussions see Refs. [44, 61, 62, 75–78].Figure 6 shows a sequence of phase-space slice plotsof regular orbits in the vicinity of the central fixed pointfor the parameter sets (A), (C), and (F), see Fig. 4. InFig. 6(a) for parameter set (A), i.e. K = 0 .
31, one isin the stable regime and quite far away from the Kreincollision. The EE fixed point (red sphere) is surroundedby regular tori shown as grey curves, which form pairsclosed loops on either side of the fixed point. The gen-eral arrangement of the tori is governed by the two(Lyapunov) families of -tori which are attached to theEE fixed point and shown in yellow and orange, respec-tively. Due to the symmetries of the map, both familieslie in the q – q plane. Thus they can be displayed in diagrams to clarify the change of the families under pa-rameter variation, see Figs. 6(d)-(f). Note that the smallgap in the yellow family in Fig. 6(a) is caused by a reso-nance, see Sec. III C. Both families of elliptic -tori aresurrounded by regular 2-tori which form pairs of ringsin the phase-space representation, depicted in graycolor. Interestingly, the regular 2-tori in the direct vicin-ity of the fixed point show a strong bending close to thefixed point. This geometry is similar to the phase spaceof the normal-form for α > I (cid:54) = 0 plane forces the tori to bendaway from the y – y plane. Furthermore, the families of tori correspond to the family of period-two periodicpoints in the normal-form.Figure 6(b) shows the situation for point (C) in pa-rameter space with K = 0 .
3. For this parameter the twoeigenvalue pairs of the linearized map at the fixed pointcoalesce at two places on the complex unit circle, seeFig. 2(b). When approaching the Krein collision param-eter, the angle between the eigenvectors of the linearizedmap decreases until the eigenvectors of the eigenvalue(a)(b)(c) p q q . . . (A) K = 0 . (d) q . . . (C) K = 0 . (e) q . . . . . . (F) K = 0 . (f) q q FIG. 6. Sequence of phase-space slice plots of regular tori represented as grey rings in the vicinity of the fixed pointshown as red spheres for elliptic-elliptic stability and as a grey sphere for complex instability. The families of tori (red,yellow, magenta) form the skeleton of the surrounding tori. The chosen parameters are (a) K = 0 .
31, (b) K = 0 . K = 0 .
285 and correspond to points (A), (C) and (F) in parameter space, see Fig. 4. The right column (d), (e), (f)depicts the families of tori, which lie in the q – q pairs become collinear. Accordingly, the families of tori are approximately parallel in the vicinity of the fixedpoint as can be seen in Fig. 6(b).Finally, Fig. 6(c) shows the situation after the Kreincollision, i.e. for K = 0 . tori detach fromthe fixed point and merge into one single family. Thiscorresponds to the normal-form behavior for α <
0, seeFig. 3(b). The regular tori close to the family of toripersist. Interestingly, orbits in the vicinity of the CUfixed point stay in its surrounding for very long timesand only eventually escape. This will be discussed inmore detail in Sec. IV.To quantify the detachment of the regular tori fromthe CU fixed point, we compute the minimal distance d tori between the complex unstable fixed point and thefamily of tori. In the normal-form description ofSec. II C the minimal distance is given by the distance be-tween the trivial fixed point at the origin and the period-two periodic point, namely by d tori = √− α/b . For the map this translates in first approximation to d tori ∝ (cid:112) K ∗ − K (24)with K ≤ K ∗ = 0 .
3. Figure 7 shows the numericallydetermined minimal distance d tori in dependence on thekicking strength K as black dots. Good agreement withthe square root behavior (24), shown as a dashed line, isfound. Further away from the Krein collision parametersmall deviations become visible.0 .
27 0 .
28 0 .
29 0 . . . . . K d tori FIG. 7. Minimal distance d tori between the central fixed pointand the family of tori in dependence on K . The distancefollows the predicted behavior ∝ (cid:112) K ∗ − K , shown as reddashed line for K ≤ K ∗ = 0 . C. Frequency space
Complementary to the representation in phase spaceone can display regular tori in frequency space, whichis particularly useful for understanding the influenceof resonances. A regular torus is characterized bytwo frequencies, one describing the motion along themajor radius of the 2-torus and one for the motionalong the minor radius. Numerically the frequencies ν , ν ∈ [0 ,
1[ for an orbit started in a phase-space point( p (0) , p (0) , q (0) , q (0)) are determined using a Fourier-transform based frequency analysis [45, 46, 79, 80]. Assignals z j ( n ) = q j ( n ) − i p j ( n ) for each degree of freedom j = 1 , q j ( n ) , p j ( n )) are the coordinatesobtained from N successive iterates of the map. In orderto distinguish regular and chaotic motion, the frequen-cies ν j of the first half of an orbit, i.e., iterates in theinterval n ∈ [0 , N/ − (cid:101) ν j of the second half, i.e., the iter-ates in the interval n ∈ [ N/ , N ]. For the motion ona regular torus, the difference of these frequency pairsshould be rather small. Thus if the maximal differencemax {| ν j − (cid:101) ν j |} is smaller than some threshold δ cut , weconsider the orbit as regular. In the following δ cut = 10 − is used. Of course, such a numerical criterion does notguarantee that the orbit eventually could become chaoticat very large times, as is also the case with other chaosindicators, see Ref. [81] for a recent overview. Using anensemble of 10 initial conditions, randomly chosen in the phase-space volume defined by p , p ∈ [ − . , .
1] and q , q ∈ [0 . , . ν , ν ) ofthe regular tori provides the two-dimensional frequencyspace representation.Figure 8 shows a sequence of such frequency space plotsfor all six parameter sets specified in Fig. 4. The fre-quencies of the EE fixed point is indicated by a large redpoint in Fig. 8(a)-(c). For the complex unstable fixed point there is only one frequency given by the angle of thecomplex eigenvalues, which is shown on the − K . As for the phase-space slice shown in Fig. 6, theorange, yellow and magenta points mark the frequenciesof the families of tori, which form the edges of thegray regions of regular tori.Resonances correspond to straight lines in frequencyspace, n ν + n ν = m, (25)with m, n , n ∈ Z and gcd( m, n , n ) = 1 and either n (cid:54) = 0 or n (cid:54) = 0. Some relevant resonance lines areshown as blue dashed lines, labeled by n : n : m . Suchresonances lead to resonance channels [46] and gaps inthe families of tori [61].The typical frequency space around an EE fixed pointis seen in Fig. 8(a)-(b) for K = 0 .
31 and K = 0 . tori are attached tothe fixed point forming a cusp and the regular tori filla region in between these families. As the eigenvaluesapproach the Krein collision parameter in Fig. 8(b), thefixed point has to approach the − ± i2 πν with ν = ν = ν . This shift ofthe frequencies of the fixed point stretches the families of tori and the top of the cusp accordingly. During thisprocess, the density of regular tori close to the − K = 0 .
3. This corre-sponds to the tangency of the families of tori so thatonly few regular tori exist in the surrounding of the fixedpoint.Figure 8(d)-(f) show the frequency space plots for thecomplex unstable case for K = 0 . K = 0 .
29, and K = 0 . tori merge in theKrein collision parameter and subsequently detach fromthe − ν , ν ) (cid:55)→ ( ν , ν − ν )to the upper branch resulting in the magenta dotted line.The transformed branch connects seamlessly to the otherbranch yielding a complete arc beginning and ending at ν ≈
0. This illustrates that both branches actually be-long to just one family of tori after the fixed pointhas turned CU. In general, such linear transformationswith determinant ± − . . .
04 (a) − : : − : : − : : − : : − : : ν (b) (c)0 .
04 0 .
06 0 . . . .
04 (d) ν ν .
04 0 .
06 0 . ν .
04 0 .
06 0 . ν FIG. 8. Frequency space for different parameters K : (a) K = 0 .
31, (b) K = 0 . K = 0 . tori. The frequency of the elliptic-elliptic (red) and complex unstable (gray) fixed point are depicted as enlargeddot. The dotted magenta curve in (d)-(f) is the unimodular transformation of the upper branch of tori. Some relevantresonance lines are shown as dashed lines. ever, the regular tori between the branches of the formercusp still exist directly after the transition as is visiblein Fig. 8(d). Only when the fixed point becomes moreunstable, the distance of the branches increases and thedensity of regular tori between them decreases until a gapemerges, see Fig. 8(f). The remaining regular orbits inFig. 8(f) are close to the family of tori. This confirmsthe observations in the phase-space slice in Fig. 6(c),where regular tori are only found in the surrounding ofthe family of tori and no regular structures are left inthe direct vicinity of the fixed point.Note that the arc like structure in the range of 0 . ≤ ν ≤ .
075 below the discussed region of regular tori, seeFig. 8(a)-(d), belongs to regular orbits in the surroundingof a periodic orbit close to the central fixed point. Al-though these orbits are not in the focus of this study theyillustrate how the complex instability of the fixed pointgradually destroys all stable structures in its vicinity.
IV. ESCAPE FROM THE CU REGION
When the EE fixed point becomes CU, this immedi-ately affects its direct surrounding as the two ellipticfamilies of tori become detached from the fixed point.Thus there are also no regular tori in its direct vicinity.Instead one has a stable and a unstable manifoldwhich lead to chaotic dynamics. However, in practiceclose to the Krein collision parameter, initial conditionsin the vicinity of the fixed point lead to orbits staying forvery long times in a confined phase-space volume. In thissection, we investigate this behavior and the underlyingescape paths in more detail. A. Escape times
To study the escape of orbits from the surroundingof the CU fixed point, we use escape time plots as inRefs. [75, 83, 84], see Fig. 9. Using a grid of initial con-ditions on a particular plane in the phase space foreach initial point the escape time n esc , required to reachsome specific exit region, is determined. Since we are in-00 . . . K = 0 . q (b)(E) K = 0 . . . . . . . K = 0 . q q . . . K = 0 . q n esc FIG. 9. Escape time plots in the q – q plane for p = 0 and p = 0 for (a) K = 0 . K = 0 .
29, (c) K = 0 . K = 0 .
27. The escape time is encoded in color, where white corresponds to those points which have not escaped within n max = 10 iterations. The fixed point is shown as red (elliptic-elliptic) or gray dot (complex unstable) and the families of tori are shown as black dots. terested in the behaviour close to the family of tori,we choose the initial points in the q – q plane throughthe fixed point with q , q ∈ [0 . , .
75] and p = 0 and p = 0. On this plane, a 2000 × q , q / ∈ [0 . , . p , p ∈ [ − . , . n max of iterations is reached.Figure 9 shows the escape time n esc encoded in colorranging from yellow for fast escaping points to black fornearly regular orbits while white points do not escape to the exit region within n max = 10 iterations (thoughthey may escape eventually). In addition, the familiesof tori are shown in black and the fixed point as redor gray dot for EE or CU stability, respectively. Theparameters for Fig. 9(a)-(c) correspond to the points (B),(E), and (F) in parameter space specified in Fig. 4. Inaddition, the families of tori are shown in black andthe fixed point as red or gray dot for EE or CU stability,respectively.As before, we focus on the structures close to the fixedpoint. For the EE case, the vicinity of the fixed point is1naturally governed by a white region which correspondsto the regular 2-tori surrounding the families of tori,compare with Fig. 9(a). Thus, even for arbitrarily largetimes these orbits do not escape. Furthermore, we seethe impact of the − K = 0 . K = 0 .
29 in Fig. 9(b) above and below thefixed point. Still, there are orbits in the direct vicinity ofthe fixed point which stay close to it for more than n max iterations. Quantitatively, the size of the white region de-pends on the threshold n max , but a larger value of n max does not affect the shown escape time plots significantly.The reason for this is that orbits in the vicinity of the CUfixed point are confined for an extremely long time whenthe parameters of the map are sufficiently close to the EEregion in Fig. 4. The more unstable the fixed point be-comes, i.e. the smaller K is, the more the two branchesof the family of tori separate and the white regiondiminishes because the regions of instability get larger.Finally, for point (F) in Fig. 4 with K = 0 .
285 all orbitsin the direct vicinity of the fixed point are able to reachthe exit region within n max iterations, see Fig. 9(c). Forthis parameter we observe that the unstable regions inthe escape plots reach the fixed point, and consequentlythe large white region is divided into two smaller ones.These two white regions correspond to the tubes of reg-ular motion in the phase-space slice representation,e.g. see Fig. 6(f), as well as the attached regular tori ofthe branches of the family of tori in frequency space,see Fig. 8(f).Figure 9(d) shows the escape time plot for K = 0 . tori moved far away from the fixed point and theunstable region in between is large. Interestingly, thisunstable region reveals a unique spiral pattern whichis attached to the fixed point. Orbits on this spiralneed at least one to two orders of magnitude more itera-tions to escape into the exit region than the neighboringones. Additionally, there is another spiral structure on asmaller scale as shown in the magnification in the inset.A closer investigation of orbits started in the darkercolored region reveals that the spiral pattern is due tothe influence of the − B. Stable and unstable manifolds
The stable and unstable manifolds associated with anunstable fixed point govern the chaotic dynamics in itssurrounding. For a complex unstable fixed point of a map the manifolds are two-dimensional invariant objectsin the phase space. Numerically the manifolds arecomputed using the parameterization method [48–50, 77,85, 86], see Appendix A for details. In the phase-space slice representation they lead to one-dimensionalcurves, see Fig. 10, where the red curve corresponds tothe unstable manifold and the blue curve to the stablemanifold.The regular 2-tori (gray loops) as well as the familiesof tori (black curves) in Fig. 10(a) are the same as inFig. 6(f). Figure 10(b) shows the geometry for a smallervalue of K = 0 .
28. The complex unstable fixed point isindicated by a gray sphere in both plots.Numerically it is found that the stable and unstablemanifolds intersect in one point. This point thereforeis a homoclinic point whose forward iterates approachthe fixed point on the stable manifold while the back-ward iterates approach the fixed point on the unstablemanifold. The existence of a transverse homoclinic pointtherefore immediately implies an infinity of such homo-clinic points. Note that generically two manifolds ina phase space will not intersect. The fact that thishappens for the manifolds of the considered fixed pointmust be due to the symmetries of the map.The geometry becomes more clearly visible for smaller K = 0 .
28 as shown in Figure 10(b). The arrangement ofthe manifolds in the phase-space slice reminds of thehomoclinic tangle in symplectic maps. In comparisonto Fig. 10(a) the excursions of the manifolds are morepronounced which corresponds to a larger chaotic regionsurrounding the complex unstable fixed point.It has to be emphasized, that even though the geome-try visually resembles the homoclinic tangle in sym-plectic maps, the iterate of any of the homoclinic inter-sections in general is not contained in the phase-spaceslice. Actually, we find numerically that the stable andunstable manifolds intersect in a line which is itselfan invariant set. Therefore, the intersection point in the phase-space slice and its iterates are only a subset ofthe intersection line. Moreover, as the manifolds areonly they cannot enclose a volume, so that there isno equivalent to the lobe structure and transport via aturnstile mechanism as in symplectic maps [87–90].2(a)(b) p q q FIG. 10. phase-space slice representation of the stable (blue) and unstable (red) manifolds of the CU fixed point togetherwith regular 2-tori (gray) and the family of tori (black) for (a) K = 0 .
285 and (b) K = 0 . .
25 0 .
26 0 .
27 0 .
28 0 .
29 0 . δ = 10 − δ = 10 − δ = 10 − K h n esc i FIG. 11. Average escape time (cid:104) n esc (cid:105) of an ensemble of 10 orbits started in U δ in dependence on K for δ = 10 − (purpledownward triangles), δ = 10 − (blue triangles), and δ = 10 − (red circles). C. Escape statistics
To investigate the chaotic transport in the vicinity ofthe CU fixed point we consider an ensemble of initialconditions in a cube U δ = [ − δ, δ ] × [0 . − δ, . δ ] , (26) with small δ . The exit region is again chosen to be p , p ∈ [ − . , .
5] and q , q / ∈ [0 . , . (cid:104) n esc (cid:105) for an ensemble of10 orbits in dependence on K for different δ = 10 − , δ = 10 − , and δ = 10 − . When approaching the Kreincollision parameter K ∗ = 0 .
3, the average escape time (cid:104) n esc (cid:105) strongly increases and for K > .
29 exceeds 10 iterations. The same is also found for the smallest escapetime (not shown). Extracting the functional dependencefrom the data turned out to inconclusive.The tail of the distribution P ( n esc ) of escape times isvery well described by an exponential, see Fig. 12. Thisprovides a hint at what mechanism could be responsiblefor such large escape times: there could be one partialbarrier (of unknown origin) for the dynamics which al-lows for a small flux towards the escape region [89]. Sucha a single partial barrier would lead to a simple exponen-tial [91] while in contrast several partial barriers wouldtypically lead to an overall power-law behavior [92–94].Note that for the small hump of (cid:104) n esc (cid:105) seen in Fig. 11around K = 0 .
284 the corresponding P ( n esc ) shows anon-exponential behavior in the tail.To quantify the escape dynamics of the ensemble, wenow consider the extent as a function of the number ofiterates. Explicitly we determine d max ( n ) = max i ≤ n {|| z ( i ) − z ∗ || | with z (0) ∈ U δ } , (27)where z ( i ) is the i -th iterate of an initial point z (0) ∈ U δ and || z i − z ∗ || is the distance to the complex unstablefixed point at z ∗ . We use 10 initial conditions in U δ with δ = 10 − . Figure 13 shows the result for five dif-ferent values of K . The expansion during the first 100iterations is similar and after about 10 iterations follows30 2 · · − − − − K = 0 . K = 0 . n esc P ( n esc ) FIG. 12. Histogram P ( n esc ) of the escape times for K =0 .
287 and K = 0 . n esc an overall exponential given by | λ | n , where λ is the eigen-value with largest absolute value. For K = 0 .
28 this isillustrated by the blue dashed curve. On a finer scale theinitial expansion happens in a step-like manner. This isdue to the spiraling motion of each orbit as illustrated inFig. 5. This motion has a different extent in the differentdirections, so that a larger distance is only obtained pe-riodically after approximately ten iterations for the firstexpansion phase. This corresponds to half the reciprocalwinding frequency of the fixed point.After the first rapid expansion phase, the maximal dis-tance shows prominent plateaus extending over several10 − − − − − − K = 0 . K = 0 . K = 0 . K = 0 . K = 0 . nd max FIG. 13. Maximal distance d max of an ensemble of10 initial conditions started in the cube U δ with δ = 10 − vs. the number of iterations n for K =0 . , . , . , . , .
29 (top to bottom, correspondingto increasing escape time). The initial expansion is well de-scribed by ∝ | λ | n , shown for K = 0 .
28 (blue dashed curve). orders of magnitude in time. These plateaus becomelonger the closer the parameter K is to K ∗ = 0 .
3, i.e.the parameter of the Krein collision. Thus for a very longtime the ensemble is effectively confined in phase space.Afterwards there is at least one trajectory which leavesthis region very quickly, as manifested by the sharp in-crease of d max .A closer look at the plateaus reveals that there is stilla rather slow increase. The occurrence of the plateauscan be explained by the alternating spiraling in and outof the dynamics already observed in Refs. [6, 11, 36, 38]:An orbit initially started near the complex unstable fixedpoint moves away from it on a spiral along the unstablemanifold until it reaches a maximal distance to the fixedpoint. This behavior corresponds to the first expansionphase up to approximately 100 iterations. Subsequently,the orbit spirals in again and gets very close to the fixedpoint with some minimal distance. When spiraling outagain, this can lead to a slightly increased maximal dis-tance. This process of inward and outward spiraling re-peats many times before the orbit escapes quickly. Note,that this sequence of outward and inward spiraling onlyholds for parameters which are near the elliptic-ellipticregion in the parameter plot in Fig. 4, i.e. if K is suffi-ciently close to K ∗ = 0 .
3. Further away from the Kreincollision parameter the extent of the plateau of d max be-comes very short or even non-existent, see Fig. 13 for K = 0 . (cid:101) d max ( n ) = max {|| z ( n ) − z ∗ || | with z (0) ∈ U δ } , (28)see Fig. 14. Initially one has the overall exponential in-crease which is superimposed by small oscillation causedby the spiraling motion. This occurs until the ensem-ble has expanded until the homoclinic intersection, whichcorresponds to the beginning of the plateaus in Fig. 13.Afterwards, there is a prominent dip around n = 200 i.e.the extent of the ensemble has become quite small againand most of the points are located in a small surroundingof the complex unstable fixed point. These minima con-verge to the plateau for growing n such that the seconddip is already barely visible. This effect is due to the in-ward and outward spiraling behavior of each individualorbit. The inset of Fig. 14 shows (cid:101) d max of one single orbit.The position of the first minimum after the expansion ofone single orbit matches roughly the first minimum in theplateau of the ensemble. This expansion and contractionof the ensemble repeats approximately periodically un-til some loss of correlations sets in and the dips of d max become less and less prominent. Note that such kindof dynamics is also found for symplectic maps for thedynamics after a period-doubling bifurcation and also for symplectic maps with an II fixed-point. A more de-tailed investigation and comparison of these cases wouldbe very interesting and is left for future studies.410 − − − − − − n e d max . . FIG. 14. Maximal extent (cid:101) d max of an ensemble of 10 initialconditions started in the cube U δ with δ = 10 − vs. thenumber of iterations n . The inset shows the maximal extentof a single exemplary orbit up to the first 1000 iterations. D. Escape dynamics
The temporal dependence of the extent of the iter-ates of the ensemble allows for quantifying the long-timeconfinement within the chaotic region surrounding thecomplex unstable fixed point. Still, the key question is,what is responsible for this long-time confinement andwhat is the escape mechanism? In particular, referringto the normal form description, there could be either anescape within the I = 0 plane or across different planeswith I (cid:54) = 0. Escape within I = 0 would be similar to thecase of the period-doubling bifurcation in symplecticmaps, where just after the fixed point has become un-stable there are usually still invariant curves so that anescape of orbits is only possible when being further awayfrom the bifurcation in parameter space. In contrast,the escape across different planes with I (cid:54) = 0 would bea genuinely effect. In principle there could also be acompetition between these two escape routes and whichof them is relevant could depend both on parameters andconsidered time-scales.As a measure of the invariant I of the normal form fora symplectic map we make use of the quadratic invariantof the linearized map. With Eq. (6) we get Q = − p + p − q ( K − K ) − q ( K − K )+ p q ( K − K ) + p q ( K − K )+ K ( p q − p q + 2 q q ) . (29)By use of a suitable coordinate transformation Eq. (29)degenerates for the Krein collision parameter into twoplanes, namely the p = − p and the q = q plane [33]. − −
50 0 50 − . . e nQ q q p q q (a)(b) 3711151923 j FIG. 15. (a) Shown are the q , q -coordinates and thequadratic invariant of the linearization Q of 1000 orbitsstarted with random initial conditions in U δ with δ = 10 − and K = 0 .
288 over (cid:101) n = n − n esc . The escape criterionis the same as in the previous experiments and marked as ared dashed line while in the last plot the blue dashed linerepresents Q . (b) phase-space slice representation of seg-ments of a single exemplary orbit for K = 0 . j = 3 , , , ,
19, and 23, see the text for further expla-nation. The unstable manifold is shown as a red curve andthe Q = Q These two planes geometrically correspond to the repre-sentation of the I = 0 plane for the hyperplanes x = 0 oralternatively y = 0 in the normal form description, seesection II C. Hence, the quadratic invariant at the fixedpoint is Q ( z ∗ ) = 0.However, away from the Krein collision the two planesare not degenerate anymore. Therefore Q does not re-semble the I = 0 plane and we get Q = Q ( z ∗ ) = − K + K K, (30)which is not zero in general. Still it turns out, that Q − Q is a well suited quantity to approximate the invariant I q and q coordi-nates as well as the quadratic invariant Q as function of (cid:101) n = n − n esc , i.e. for a few iterations before and afterthe escape of an orbit. The initial conditions of the en-semble with 1000 orbits are started in U δ with δ = 10 − and the kicking strength is K = 0 . (cid:101) n and fulfill the escape criterionfor positive (cid:101) n , as indicated by the red dashed horizontallines. The spread of the distances of q and q aroundthe fixed point, i.e. the width of the distribution of dis-tances around 0.5, is slightly increasing towards (cid:101) n = 0.Even though this trend is visible in both coordinates, theescape condition is reached first by the q coordinate.In order to understand the escape mechanism in termsof the phase space geometry, we compare the escape pathin the phase-space slice with the geometry of thenormal-form. The arrangement of regular tori and thefamily of tori, see Fig. 6, suggest that the Q = Q plane is a good approximation to the I = 0 plane, com-pare to the gray plane in Fig. 15(b). Therefore, Q − Q provides an approximate measure of how far a point of anorbit is away from the I = 0 plane, see Fig. 15(a). As forthe single coordinates, Q shows an overall increase and isspread more widely as (cid:101) n approaches 0. However, about7 iterations before (cid:101) n = 0 the distribution of Q splits intotwo separate parts, away from 0.In order to determine if the ensemble escapes throughthese two separated escape paths or interchanges betweenthose two, we split the ensemble in two subsets by either Q ( (cid:101) n = 0) > Q or Q ( (cid:101) n = 0) < Q and determine theirmean and variance. Figure 16(a) shows the average asdots and their standard deviation as error bars of the Q ( (cid:101) n = 0) > Q and the Q ( (cid:101) n = 0) < Q subset in blueand red color, respectively. The ensemble clearly sepa-rates in these two sets and fluctuates around Q markedas a black dashed line. Once the escape criterion is ful-filled, either Q > Q or Q < Q and initially no fur-ther change in sign occurs. This behavior translates toescape either across I >
I < I is onlypossible because the normal-form geometry provided byEq. (18) is broken.Figure 16(b) shows the time evolution of the varianceof both sets ranging from 2000 iterations before the es-cape up to the escape. We observe the same type of in-crease of the variance for both subsets towards the escapeat (cid:101) n = 0. Understanding the behaviour of the variancequantitatively is an interesting future task.By following one single orbit we can also get an intu-ition of how the orbit crosses the different I (cid:54) = 0 planes,see Fig. 15(b). Here, we consider a single orbit withinitial condition ( p , p , q , q ) = (0 , , . µ, . µ )with µ = 10 − for K = 0 . − −
50 0 − . . . e n h Q i a) − − − − e n var( Q )b) FIG. 16. The mean width standard deviation a) and thevariance b) of the quadratic invariants Q of the ensemble inFig. 15 are shown. The blue set corresponds to the set oforbits with Q ( (cid:101) n = 0) > Q and the inverse to the red datapoints. The black dashed line in a) represents Q mentation of the map. For this orbit we consider succes-sive segments [ j · , ( j +1) · ε = 10 − are determined. A selectionin the surrounding of the complex unstable fixed pointis shown in Fig. 15(b) together with the phase-spaceslice of the unstable manifold as a red curve. This plotshows that the iterates of the initial point are approxi-mately restricted around lines in the phase-spaceslice. These lines follow the unstable manifold and eachof the successive segments appears to lie on a slightlybent surface, similar to the I (cid:54) = 0 planes of the normalform, compare with Fig. 3. This suggests that an escap-ing orbit is following the unstable manifold which givesrise to transport through the I (cid:54) = 0 planes. Note that theslice segments for j = 7 ,
19 are located at the excursion ofthe manifold farther away from the fixed point and there-fore do only appear at the edge of the magnification. Ingeneral, the motion along the unstable manifold explainsalso the repetitive expanding and contracting behavior ofthe orbits.6
V. SUMMARY AND OUTLOOK
In this paper the transition of a fixed point withelliptic-elliptic dynamics to complex-unstable dynamicsunder parameter variation is investigated for a sym-plectic map. Using phase-space slices we visualizeregular dynamics in the vicinity of the fixed point. Whilein the elliptic-elliptic case there exist two families of tori which are attached to the fixed point and are sur-rounded by regular 2-tori, these families merge into onesingle family and split off the fixed point. Moreover, thegeometry of regular orbits close to the fixed point in the phase-space slice lie on surfaces as predicted by thenormal form description, see Fig. 3. The phase-space rep-resentation is complemented by a frequency analysis ofregular tori, see Fig. 8. Before the transition to complexinstability the two families of tori are attached to thefixed point forming a cusp-like region which encloses theregular tori. The fixed point becomes complex unstableunder parameter variation when reaching the − tori split off thefixed point. Applying a unimodular transformation clar-ifies that these apparently two families of tori actuallyform a single arc in frequency space.Once the fixed point has become complex unstablenearby orbits may eventually escape. However, it turnsout that shortly after the transition orbits are confinedto a particular phase-space region for very long times.This region can be visualized using escape time plots,see Fig. 9. The extent is governed by the stable and un-stable invariant manifolds of the complex unstable fixedpoint. In the phase space they lead to a geometrywhich is visually similar to that of the well-known homo-clinic tangle for symplectic maps.To quantify these observations we consider the escapestatistics for an ensemble of 10 orbits, started in thevicinity of the fixed point in dependence on the distanceto the bifurcation point, i.e. by varying the parameter K . The average escape time strongly increases whenapproaching the bifurcation point. Measuring the max-imal distance of all orbits of the ensemble to the fixedpoint over the number of iterations, reveals three differ-ent phases of the dynamics, see Fig. 13. Initially, for thefirst approximately 100-200 iterations, the distance in-creases exponentially, followed by a extended plateaus inthe second phase. These plateaus correspond to the long-time confinement and extend over longer times the closerthe parameter is to the Krein bifurcation. A closer lookat the plateaus shows that there is a very slow increaseas function of time. The plateaus are due to the inwardand outward spiraling dynamics of the ensemble whichfollows the unstable invariant manifold. Thus, the slopecorresponds to a slowly growing extent of individual or-bits, see Fig. 14. Eventually, in the last phase one orbitof the ensemble will escape after a critical time and themaximal distance of the ensemble quickly reaches approx-imately 1. If the fixed point is very unstable, the plateauis very short or even not existent. Comparing the q , q coordinates and the quadraticinvariant Q of the ensemble for the transition from phasetwo to three allows for determining the main escape pathsclose to the bifurcation, see Fig. 15. This provides evi-dence that long confined orbits escape across either I >
I < I for the specific map using a numerical normalform analysis. This would allow for accurately quanti-fying the transport across the approximately invariantplanes and to investigate whether the escape can be de-scribed by a diffusive process. ACKNOWLEDGMENTS
We are grateful for discussions with Markus Firm-bach, Franziska H¨ubner, Roland Ketzmerick, and HarisSkokos. Robert MacKay kindly provided us with acopy of Ref. [40]. Furthermore, we acknowledge supportby the Deutsche Forschungsgemeinschaft under grantKE 537/6–1.All visualizations were created using Mayavi [95].
Appendix A: Computing stable and unstablemanifolds
There are various methods to determine the invariantmanifolds associated with an unstable fixed point, see e.g.Refs. [96–101]. Here we use the parametrization method,which was introduced in Refs. [48–50] and has been usedfor example in Refs. [77, 85, 86].The parameterization method takes advantage of theHartman-Grobman theorem which for symplectic mapsstates that the linearization of a fixed point or periodicorbit is conjugate to its local stable and unstable in-variant manifolds W s , uloc , if the eigenvalues have an ab-solute value different from one, i.e., if they are unstable.The key point of the parameterization method is to findsmooth vector-valued functions F s and F u which param-eterize the stable and unstable invariant manifolds. Inorder to do so, F s , u have to obey on the one hand thelinear conditions F s,u ( ) = z ∗ and (A1) ∂ F s,u ( θ ) ∂θ j = ξ j for 1 ≤ i ≤ n s,u (A2)with θ = ( θ , . . . , θ n s,u ) ∈ C n s,u and ξ j ∈ C n s,u be-ing the associated eigenvector to the n s stable and n u unstable eigenvalues λ j .On the other hand, F s,u must satisfy the conjugacyequation M ◦ F s,u ( θ ) = F s,u ( λ θ , . . . , λ n s,u θ n s,u ) (A3)7in order to take the non-linearity of the map into account.For a complex unstable fixed point of a symeplecticmap one finds n s = n u = 2. Therefore, we expand F s,u into the power series F s,u ( θ , θ ) = p ( θ , θ ) p ( θ , θ ) q ( θ , θ ) q ( θ , θ ) = ∞ (cid:88) i =0 ∞ (cid:88) j =0 f ij θ i θ j , (A4)with vector-valued coefficients f ij ∈ C .For the considered map (18), the non-linear terms consist of sine functions with various input arguments,namely three different sums of phase-space coordinates.We approximate these sine functions by their Tay-lor series representation. Advantageously, the coeffi-cients of this series can be easily computed by an auto-differentiation algorithm which is based on Refs. [102,103]. Using the series representation of the sine termsof the map M and combining (A4) and the conjugacyequation (A3) leads to a homological equation which canbe solved iteratively for the coefficients f ij up to a givenorder ( m, n ). The corresponding initial value problem issolved by the linear conditions (A1) and (A2). [1] R. Broucke, Stability of periodic orbits in the elliptic,restricted three-body problem. , AIAA Journal , 1003(1969).[2] J. D. Hadjidemetriou, The stability of periodic orbits inthe three-body problem , Celestial Mech. , 255 (1975).[3] B. Sicardy, Stability of the triangular Lagrange pointsbeyond Gascheau’s value , Celest. Mech. Dyn. Astron. , 145 (2010).[4] P. Magnenat,
Numerical study of periodic orbit proper-ties in a dynamical system with three degrees of freedom ,Celestial Mech. , 319 (1982).[5] L. Martinet and D. Pfenniger, Complex instabilityaround the rotation axis of stellar systems. I. Galacticpotentials , Astron. & Astrophys. , 81 (1987).[6] D. C. Heggie,
Bifurcation at complex instability , Celes-tial Mech. , 357 (1985).[7] G. Contopoulos and P. Magnenat, Simple three-dimensional periodic orbits in a galactic-type potential ,Celestial Mech. , 387 (1985).[8] P. A. Patsis and L. Zachilas, Complex instability of sim-ple periodic orbits in a realistic two-component galacticpotential , Astron. & Astrophys. , 37 (1990).[9] M. Oll´e, J. R. Pacha, and J. Villanueva,
Motion closeto the hopf bifurcation of the vertical family of periodicorbits of l4 , Celest. Mech. Dyn. Astron. , 87 (2004).[10] H. Hanßmann and J.-C. van der Meer, Algebraic meth-ods for determining Hamiltonian Hopf bifurcations inthree-degree-of-freedom systems , J. Dyn. Diff. Equat. ,455 (2005).[11] M. Katsanikas, P. A. Patsis, and G. Contopoulos, Thestructure and evolution of confined tori near a Hamil-tonian Hopf bifurcation , Int. J. Bifurcation Chaos ,2321 (2011).[12] P. A. Patsis and M. Katsanikas, The phase spaceof boxy-peanut and X-shaped bulges in galaxies — II.The relation between face-on and edge-on boxiness ,Mon. Not. R. Astron. Soc. , 3546 (2014).[13] K. Efstathiou, R. H. Cushman, and D. A. Sadovski´ı,
Hamiltonian Hopf bifurcation of the hydrogen atom incrossed fields , Physica D , 250 (2004).[14] A. Lahiri and M. S. Roy,
The Hamiltonian Hopf bifur-cation , International Journal of Non-Linear Mechanics , 787 (2001).[15] M. Oll´e and J. R. Pacha, Hopf bifurcation for the hy-drogen atom in a circularly polarized microwave field ,Commun. Nonlinear Sci. Numer. Simulat. , 27 (2018). [16] S. C. Farantos and M. Founargiotakis, Wave packetdynamics and phase space structure of HCN molecule ,Chemical Physics , 345 (1990).[17] G. D´ıaz, J. Egea, S. Ferrer, J. C. van der Meer, andJ. A. Vera,
Relative equilibria and bifurcations in thegeneralized van der Waals 4D oscillator , Physica D ,1610 (2010).[18] J. E. Howard, A. J. Lichtenberg, M. A. Lieberman, andR. H. Cohen,
Four-dimensional mapping model for two-frequency electron cyclotron resonance heating , Phys-ica D , 259 (1986).[19] J.-C. van der Meer, The Hamiltonian Hopf Bifurcation ,number 1160 in Lecture Notes in Mathematics, SpringerBerlin Heidelberg (1985).[20] J. D. Crawford,
Introduction to bifurcation theory ,Rev. Mod. Phys. , 991 (1991).[21] H. Papadaki, G. Contopoulos, and C. Polymilis, Com-plex instability , in A. E. Roy and B. A. Steves (editors)“From Newton to Chaos”, 485, Plenum Press, New York(1995).[22] P. D. McSwiggen and K. R. Meyer,
The evolution ofinvariant manifolds in Hamiltonian-Hopf bifurcations ,J. Diff. Eqs. , 538 (2003).[23] M. Oll´e, J. R. Pacha, and J. Villanueva,
Dynamics andbifurcation near the transition from stability to complexinstability , in J. Delgado, E. A. Lacomba, J. Llibre, andE. P´erez-Chavela (editors) “New Advances in CelestialMechanics and Hamiltonian Systems”, 185, Springer US(2004).[24] M. Oll´e, J. R. Pacha, and J. Villanueva,
Dynamics closeto a non semi-simple 1:-1 resonant periodic orbit , Dis-crete Contin. Dyn. Sys. Ser. B , 799 (2005).[25] E. Fontich, C. Sim´o, and A. Vieiro, Splitting of the sepa-ratrices after a Hamiltonian–Hopf bifurcation under pe-riodic forcing , Nonlinearity , 1440 (2019).[26] G. Wen, Criterion to identify Hopf bifurcations in mapsof arbitrary dimension , Phys. Rev. E , 026201 (2005).[27] H. W. Broer, H. Hanßmann, and J. Hoo, The quasi-periodic Hamiltonian Hopf bifurcation , Nonlinearity ,417 (2007).[28] M. Oll´e, J. R. Pacha, and J. Villanueva, Kolmogorov–Arnold–Moser aspects of the periodic hamiltonian Hopfbifurcation , Nonlinearity , 1759 (2008).[29] R. Vitolo, H. Broer, and C. Sim´o, Quasi-periodic bifur-cations of invariant circles in low-dimensional dissipa- tive dynamical systems , Regul. Chaotic Dyn. , 154(2011).[30] G. Contopoulos, S. C. Farantos, H. Papadaki, andC. Polymilis, Complex unstable periodic orbits andtheir manifestation in classical and quantum dynamics ,Phys. Rev. E , 4399 (1994).[31] J. E. Howard and R. S. MacKay, Linear stability of sym-plectic maps , J. Math. Phys. , 1036 (1987).[32] Ch. Skokos, On the stability of periodic orbits of high di-mensional autonomous Hamiltonian systems , Physica D , 155 (2001).[33] D. Pfenniger,
Numerical study of complex instability. I.Mappings , Astron. & Astrophys. , 97 (1985).[34] D. Pfenniger,
Numerical study of complex instability.II. barred galaxy bulges , Astron. & Astrophys. , 112(1985).[35] C. Froeschle,
On the number of isolating integrals insystems with three degrees of freedom , Astrophys. SpaceSci. , 110 (1971).[36] `A. Jorba and M. Oll´e, Invariant curves nearHamiltonian-Hopf bifurcations of four-dimensional sym-plectic maps , Nonlinearity , 691 (2004).[37] L. Zachilas, M. Katsanikas, and P. A. Patsis, The struc-ture of phase space close to fixed points in a 4D symplec-tic map , Int. J. Bifurcation Chaos , 1330023 (2013).[38] N. Delis and G. Contopoulos, Analytical and nu-merical manifolds in a symplectic 4-D map , Ce-lest. Mech. Dyn. Astron. , 313337 (2016).[39] T. K. Roy and A. Lahiri,
Reversible Hopf bifurcation infour-dimensional maps , Phys. Rev. A , 4937 (1991).[40] T. J. Bridges, Cushman, R. H., and R. S. MacKay, Dy-namics near an irrational collision of eigenvalues forsymplectic maps , Fields Institute Communications , 61(1995).[41] E. Fontich, C. Sim´o, and A. Vieiro, The discreteHamiltonian–Hopf bifurcation for 4D symplectic maps ,in M. Corbera, J. M. Cors, J. Llibre, and A. Ko-robeinikov (editors) “Extended Abstracts Spring 2014”,Number 4 in Trends in Mathematics, 77, Springer In-ternational Publishing Switzerland (2015).[42] A. Lahiri, A. Bhowal, and T. K. Roy,
Resonant colli-sions in four-dimensional reversible maps: A descrip-tion of scenarios , Physica D , 95 (1998).[43] A. Bhowal, T. K. Roy, and A. Lahiri,
Hopf bifurcation infour-dimensional reversible maps and renormalisationequations , Phys. Lett. A , 9 (1993).[44] M. Richter, S. Lange, A. B¨acker, and R. Ketzmerick,
Visualization and comparison of classical structures andquantum states of four-dimensional maps , Phys. Rev. E , 022902 (2014).[45] J. Laskar, The chaotic motion of the solar system: A nu-merical estimate of the size of the chaotic zones , Icarus , 266 (1990).[46] J. Laskar, Frequency analysis for multi-dimensional sys-tems. Global dynamics and diffusion , Physica D , 257(1993).[47] J. Laskar, Frequency analysis of a dynamical system ,Celest. Mech. Dyn. Astron. , 191 (1993).[48] X. Cabr´e, E. Fontich, and R. de la Llave, The parameter-ization method for invariant manifolds I: Manifolds as-sociated to non-resonant subspaces , Indiana Univ. Math.J. , 283 (2003). [49] X. Cabr´e, E. Fontich, and R. de la Llave, The parame-terization method for invariant manifolds II: Regularitywith respect to parameters , Indiana Univ. Math. J. ,329 (2003).[50] X. Cabr´e, E. Fontich, and R. de la Llave, The parameter-ization method for invariant manifolds III , J. Diff. Eqs. , 444 (2005).[51] J. Moser,
New aspects in the theory of stability of Hamil-tonian systems , Comm. Pure Appl. Math. , 81 (1958).[52] V. I. Arnold and A. Avez, Ergodic Problems of ClassicalMechanics , Benjamin, NewYork (1968).[53] M. G. Kre˘ın,
Topics in Differential and Integral Equa-tions and Operator Theory , number 7 in Operator The-ory: Advances and Applications (edited by I. Gohberg),Birkh¨auser, Basel (1983).[54] T. J. Bridges and R. H. Cushman,
Unipotent normalforms for symplectic maps , Physica D , 211 (1993).[55] T. J. Bridges and J. E. Furter, Singularity Theoryand Equivariant Symplectic Maps , number 1558 in Lec-ture Notes in Mathematics, Springer Berlin Heidelberg(1993).[56] K. Meyer and D. C. Offin,
Introduction to HamiltonianDynamical Systems and the N-Body Problem , SpringerInternational Publishing, Cham, third edition (2017).[57] S. M. Graff,
On the conservation of hyperbolic invarianttori for Hamiltonian systems , J. Diff. Eqs. , 1 (1974).[58] E. Zehnder, Generalized implicit function theoremswith applications to some small divisor problems, II ,Comm. Pure Appl. Math. , 49 (1976).[59] `A. Jorba and J. Villanueva, On the normal behaviour ofpartially elliptic lower-dimensional tori of Hamiltoniansystems , Nonlinearity , 783 (1997).[60] `A. Jorba and J. Villanueva, The fine geometry of theCantor families of invariant tori in Hamiltonian sys-tems , in C. Casacuberta, R. Mir´o-Roig, J. Verdera, andS. Xamb´o-Descamps (editors) “European Congress ofMathematics”, volume 202 of
Progress in Mathematics ,557, Birkh¨auser Basel (2001).[61] S. Lange, M. Richter, F. Onken, A. B¨acker, and R. Ket-zmerick,
Global structure of regular tori in a generic 4Dsymplectic map , Chaos , 024409 (2014).[62] F. Onken, S. Lange, R. Ketzmerick, and A. B¨acker, Bi-furcations of families of 1D-tori in 4D symplectic maps ,Chaos , 063124 (2016).[63] G. Contopoulos, Qualitative changes in 3-dimensionaldynamical systems , Astron. & Astrophys. , 244(1986).[64] G. Contopoulos and A. Giorgilli,
Bifurcations and com-plex instability in a 4-dimensional symplectic mapping ,Meccanica , 19 (1988).[65] T. M. Cherry, On periodic solutions of Hamiltonian sys-tems of differential equations , Phil. Trans. R. Soc. A , 137 (1928).[66] K. R. Meyer,
Generic bifurcation of periodic points ,Trans. Amer. Math. Soc. , 95 (1970).[67] J. M. Greene, R. S. MacKay, F. Vivaldi, and M. J.Feigenbaum,
Universal behaviour in families of area-preserving maps , Physica D , 468 (1981).[68] R. S. MacKay, Renormalisation in area–preservingmaps , number 6 in Advanced Series in Nonlinear Dy-namics, World Scientific, Singapure (1993). [69] K. Meyer, G. Hall, and D. Offin, Introduction to Hamil-tonian Dynamical Systems and the N-Body Problem ,Springer-Verlag, New York (2009).[70] C. Froeschl´e,
Numerical study of dynamical systemswith three degrees of freedom. I. Graphical displays offour-dimensional sections , Astron. & Astrophys. , 115(1970).[71] C. Froeschl´e, Numerical study of a four-dimensionalmapping , Astron. & Astrophys. , 172 (1972).[72] `A. Haro, Center and center-(un)stable manifolds ofelliptic-hyperbolic fixed points of 4d-symplectic maps.an example: The froeschl´e map , in C. Sim´o (editor)“Hamiltonian Systems with Three or More Degrees ofFreedom”, volume 533 of
NATO ASI Series: C - Math-ematical and Physical Sciences , 403, Kluwer AcademicPublishers, Dordrecht (1999).[73] M. Guzzo, E. Lega, and C. Froeschl´e,
On the numericaldetection of the effective stability of chaotic motions inquasi-integrable systems , Physica D , 1 (2002).[74] A. Celletti, C. Falcolini, and U. Locatelli,
On the break-down threshold of invariant tori in four dimensionalmaps , Regular & Chaotic Dynamics , 227 (2004).[75] A. B¨acker and J. D. Meiss, Elliptic bubbles inMoser’s 4D quadratic map: The quadfurcation , SIAMJ. Appl. Dyn. Syst. , 442 (2020).[76] S. Lange, A. B¨acker, and R. Ketzmerick, What isthe mechanism of power-law distributed Poincar´e recur-rences in higher-dimensional systems? , EPL , 30002(2016).[77] S. Anastassiou, T. Bountis, and A. B¨acker,
Homo-clinic points of 2D and 4D maps via the parametrizationmethod , Nonlinearity , 3799 (2017).[78] M. Firmbach, S. Lange, R. Ketzmerick, and A. B¨acker, Three-dimensional billiards: Visualization of regu-lar structures and trapping of chaotic trajectories ,Phys. Rev. E , 022214 (2018).[79] C. C. Martens, M. J. Davis, and G. S. Ezra, Local fre-quency analysis of chaotic motion in multidimensionalsystems: Energy transport and bottlenecks in planarOCS , Chem. Phys. Lett. , 519 (1987).[80] R. Bartolini, A. Bazzani, M. Giovannozzi, W. Scandale,and E. Todesco,
Tune evaluation in simulations and ex-periments , Part. Accel. , 147 (1996).[81] Ch. Skokos, G. A. Gottwald, and J. Laskar (edi-tors) Chaos Detection and Predictability , volume 915of
Lecture Notes in Physics , Springer Berlin Heidelberg(2016).[82] M. Born,
The Mechanics of the Atom , G. Bell & Sons,Ltd., London (1927).[83] P. Gaspard and S. A. Rice,
Hamiltonian mapping modelsof molecular fragmentation , J. Phys. Chem. , 6947(1989).[84] T. Schilling,
4D phase space and escape in van derWaals molecules , Masterthesis, Technische Universit¨atDresden, Fachrichtung Physik (2017).[85] `A. Haro, M. Canadell, J.-L. Figueras, A. Luque, and J.-M. Mondelo,
The Parameterization method for invari- ant manifolds: from rigorous results to effective compu-tations , volume 195 of
Applied Mathematical Sciences ,Springer International Publishing (2016).[86] J. Gonzalez and J. D. Mireles James,
High-order pa-rameterization of stable/unstable manifolds for long pe-riodic orbits of maps , SIAM J. Appl. Dyn. Syst. ,1748 (2017).[87] R. S. MacKay, J. D. Meiss, and I. C. Percival, Transportin Hamiltonian systems , Physica D , 55 (1984).[88] V. Rom-Kedar and S. Wiggins, Transport in two-dimensional maps , Arch. Rational Mech. Anal. , 239(1990).[89] J. D. Meiss,
Symplectic maps, variational principles,and transport , Rev. Mod. Phys. , 795 (1992).[90] J. D. Meiss, Thirty years of turnstiles and transport ,Chaos , 097602 (2015).[91] W. Bauer and G. F. Bertsch, Decay of ordered andchaotic systems , Phys. Rev. Lett. , 2213 (1990).[92] M. Ding, T. Bountis, and E. Ott, Algebraic escape inhigher dimensional Hamiltonian systems , Phys. Lett. A , 395 (1990).[93] B. V. Chirikov and V. V. Vecheslavov,
Theoryof fast Arnold diffusion in many-frequency systems ,J. Stat. Phys. , 243 (1993).[94] D. L. Shepelyansky, Poincar´e recurrences in Hamilto-nian systems with a few degrees of freedom , Phys. Rev. E , 055202(R) (2010).[95] P. Ramachandran and G. Varoquaux, Mayavi: 3D vi-sualization of scientific data , Comput. Sci. Eng. , 40(2011).[96] Y. You, E. J. Kostelich, and J. A. Yorke, Calculatingstable and unstable manifolds , Int. J. Bifurcation Chaos , 605 (1991).[97] D. Hobson, An efficient method for computing invariantmanifolds of planar maps , J. Comput. Phys. , 14(1993).[98] E. J. Kostelich, J. Yorke, and Z. You,
Plotting sta-ble manifolds: error estimates and noninvertible maps ,Physica D , 210 (1996).[99] R. H. Goodman and J. K. Wr´obel, High-order bisec-tion method for computing invariant manifolds of two-dimensional maps , Int. J. Bifurcation Chaos , 2017(2011).[100] J. K. Wr´obel and R. H. Goodman, High-order adap-tive method for computing two-dimensional invariantmanifolds of three-dimensional maps , Commun. Non-linear Sci. Numer. Simulat. , 1734 (2013).[101] C. Efthymiopoulos, G. Contopoulos, and M. Kat-sanikas, Analytical invariant manifolds near un-stable points and the structure of chaos , Ce-lest. Mech. Dyn. Astron. , 331 (2014).[102] R. D. Neidinger,