Regular Regimes of the Three Body Harmonic System

Regular Regimes of the Three Bo dy Harmonic System Ori Sap orta Katz Dep artment of Applie d Mathematics, Weizmann Institute of Scienc e, R ehovot 76100, Isr e al Efi Efrati ∗ Dep artment of Physics of Complex Systems, Weizmann Institute of Scienc e, R ehovot 76100, Isr ael (Dated: December 18, 2019) The symmetric harmonic three-mass system with finite rest lengths, despite its apparent sim- plicit y , displa ys a wide array of interesting dynamics for different energy v alues. A t low energy the system sho ws regular b eha vior that pro duces a deformation-induced rotation with a constant a veraged angular v elo city . As the energy is increased this b eha vior makes w ay to a chaotic regime with rotational b ehavior statistically resembling L ´ evy walks and random walks. A t high enough energies, where the rest lengths b ecome negligible, the chaotic signature v anishes and the system returns to regularity , with a single dominan t frequency . The transition to and from chaos, as well as the anomalous pow er law statistics measured for the angular displacement of the harmonic three mass system are largely gov erned b y the structure of regular solutions of this mixed Hamiltonian system. Th us a deeper understating of the system’s irregular b eha vior requires mapping out its regular solutions. In this w ork we pro vide a comprehensive analysis of the system’s regular regimes of motion, using p erturbativ e metho ds to derive analytical expressions of the system as almost-in tegrable in its lo w- and high-energy extremes. The compatibility of this description with the full system is shown n umerically . In the lo w-energy regime, the Birkhoff normal form metho d is utilized to circumv ent the lo w-order 1:1 resonance of the system, and the conditions for Kolmogoro v-Arnold-Moser theory are sho wn to hold. The in tegrable appro ximations provide the back-bone structure around which the b ehavior of the full non-linear system is organized, and pro vide a path wa y to understanding the origin of the p ow er-law statistics measured in the system. I. INTR ODUCTION Recen tly , the harmonic three-mass system with finite rest lengths (Fig. 1), was studied and its statistical b e- ha vior analyzed [1]. This deceptively simple system was sho wn to display a ric h v ariet y of dynamics due to geo- metric non-linearities induced b y the finite rest lengths of the springs, which render the system dynamically mixed. F or different energies and initial conditions, the system exhibits constan t deformation-induced rotation with zero angular momen tum and random w alk of the orien tation angle, among other phenomena. P erhaps the most sur- prising dynamical feature exhibited by the system is the L ´ evy-w alk regime: for a con tinuous range of energies, the orien tation of the system as a rotating triangle p er- forms b outs of constant av erage velocity , switc hing di- rections with a p o wer-la w distribution, fitting the L´ evy w alk mo del [2 – 4]. The anomalous exp onent attributed to this dynamics seems to in terp olate smo othly b et ween the v alue of 2, signifying coherent ballistic b ehavior, and 1, signifying regular random w alk statistics. In lo w-dimensional systems, d ≤ 2, the emergence of p ow er-law statistics is well-understoo d, attributed to the breakdo wn of Kolmogoro v-Arnold-Moser (KAM) tori creating partial transp ort barriers which can b e crossed b y chaotic tra jectories at a slow rate in a phenomenon commonly referred to as sticking or trapping [5 – 10]. ∗ efi.efrati@weizmann.ac.il Ho wev er, despite the robustness of this phenomenon [4, 11 – 14], a general framework for the origin of p ow er- la w statistics in high-dimensional mixed Hamiltonian sys- tems con tin ues to elude current understanding [5, 15–17]. In the three-b o dy harmonic system, the coherent b outs creating the p ow er-la w statistics strongly resemble their lo wer-energy regular counterparts, indicating a partial trapping of chaotic tra jectories around regular islands for finite times. A quantitativ e analysis of the dynamical mec hanism b ehind this phenomenon would require a deep understanding of the regular b ehavior of the non-linear system. In this work we seek to identify and characterize the regular solutions of the harmonic three b ody system, complemen ting the work in [1], in the extreme low- and high-energy regimes. By using a p erturbative approach w e find integrable appro ximations of the Hamiltonian and characterize their solutions. W e sho w how presenting the dynamics of the full system in the phase space v ari- ables induced by the integrable approximations leads to a simplified picture that allows a clearer interpretation of the c haotic dynamics. Section II presents the system and its interesting dynamics, and provides a brief summary of the results of [1], as well as some extensions. Section II I deals with the low-energy regular motion, where energy confines springs to small oscillations. The Birkhoff nor- mal form metho d is employ ed in order to obtain a faith- ful description of the full system as an almost-integrable system, and the conditions for the Kolmogorov-Arnold- Moser (KAM) theory are shown to hold. In Section IV w e analyze the high-energy regular motion, where the 2 springs’ rest lengths b ecome practically negligible and the system b eha ves like the harmonic three mass system with v anishing rest lengths (which is quadratic and thus in tegrable). Section V contains a short summary and a discussion of the outlo ok of this w ork. This analysis sets the stage for a more complete under- standing of the b eha vior observed for intermediate ener- gies, as the phase space structure in the regular regimes is somewhat retained in the anomalous regimes close enough to the regular regimes, and a gradual breaking of this structure results in a contin uous transition to, and from, fully chaotic b ehavior as the energy is raised. I I. THE HARMONIC THREE MASS SYSTEM The Hamiltonian of the planar, fully symmetric three- mass system with non-zero rest lengths is H = 3 X i =1 p 2 i 2 m + X k 2 ( r ij − L ) 2 , (1) where r i = ( x i , y i ), r ij = r i − r j and r ij = | r ij | ≡ √ r ij · r ij for i, j = 1 , 2 , 3. The mass m , spring constan t k and rest length L give rise to a natural time scale τ s = p m/k , and energy scale E s = 3 2 k L 2 , the energy it takes to con tract the system to a p oin t. The parameters we use in sim ulations and in the following calculations are L = 2, k = 1 and m = 1, giving the typical time τ s = 1 and natural energy scale E s = 6. Zero-energy equilibrium is ac hieved when the distances b etw een the masses equal the rest lengths and the masses are at rest. Conserv ation of linear and angular momentum reduce the 12-dimensional phase space of the system to a 6- dimensional phase space. Energy conserv ation further reduces the dimension of the submanifold of any given tra jectory to five. As the internal motions and center-of- mass motion of the system decouple, it is straightforw ard to eliminate the four center of mass co ordinates from the Hamiltonian. The reduction of the t wo angular momen- tum degrees of freedom requires Routh reduction [18] due to the non-holonomic nature of the conserv ation law, as describ ed in section I I I. As a result, even when setting the ov erall angular momentum of the syste m to zero, the distorting triangle may exhibit deformation-induced ro- tation, appearing as a manifestation of a relev an t geomet- ric phase [19 – 21]. Thus, the orien tation of the triangle is a non-trivial, history-dep enden t v ariable of the system, and serves as a sensitive measurable for the type of dy- namics the system follows. Indeed, the dynamics of the system is incredibly rich: despite the fully harmonic interactions, the non-zero rest lengths of the springs contribute geometrical non- linearities to the system, as can b e seen algebraically in the square-ro ot term of the p otential energy . This ren- ders the system dynamically mixed, with regions of reg- ular and chaotic dynamics. Reducing our scop e to only m 1 ! ρ 1 ˆ x θ (a) 3 3 / 2 3 / 2 (b) (c) (d) m 3 m 2 ! ρ 2 φ FIG. 1. (a) The symmetric harmonic 3-mass system, with equal spring rest lengths L , equal spring constants k and equal masses m . ~ ρ 1 and ~ ρ 2 are the mass-weigh ted Jacobi co ordinates, φ is the angle b et ween them and θ is the ori- en tation v ariable of the triangle. (b), (c) and (d) are the system’s normal modes and corresp onding frequencies, com- monly known as the symmetric stretch, isometric bend and asymmetric b end, resp ectiv ely . zero angular momentum configurations, we can c harac- terize the dynamical regime by observing the orientation dynamics. In Fig. 2, the orientation of the triangle is sho wn for different regimes of motion, and the underly- ing character of the 6-dimensional dynamics, be it regu- lar, anomalous or chaotic, is apparent through this one- dimensional measurable. W e note that the numerics in Fig. 2 and throughout this work hav e b een done using the symplectic integrator provided in [22], using the sym- plectic Euler metho d [23]. A random exploration of different initial conditions sho ws that for a large p ortion of tra jectories the total energy of the system suffices to describ e the statistical qualit y of the dynamics (see Figure 3, and [1]), despite the complex structure of the mixed phase space. F or v ery low energy v alues 0 < E  E s (Fig. 2(a,b)), as w ell as for v ery high energy v alues E s ≪ E (Fig. 2(f )), the system displa ys stable quasi-p erio dic regular tra jec- tories with constant av eraged deformation-induced rota- tion rates, and a v anishing Ly apunov exp onent. In the range of energies E s / 9 . E . 2 E s (Fig. 2(d)), orienta- tion tra jectories statistically resemble uncorrelated ran- dom w alks, with a squared mean angular displacement exp onen t of 1, and a corresp onding p ositiv e Lyapuno v exp onen t. It is in the transition b etw een these three regimes that the exotic dynamics of this system is apparen t. F or energy v alues in the range E s / 15 . E . E s / 9 (Fig. 2(c)), most tra jectories exhibit a p ositiv e Lyapuno v ex- p onen t, signifying chaotic dynamics. How ever, the corre- sp onding squared mean angular displacement exp onent is anomalous, transitioning smo othly from the v alue 2, 3 (a) (b) (c) (d) (e) (f) * * * * linear trend subtracted FIG. 2. T ypical angular tra jectories of the system for v arious energies. (a) with E = 0 . 005, and (b) with E = 0 . 224, b oth exhibit regular behavior for practically infinite times. (a) has a single dominan t frequency p 3 / 2, which is the twice-degenerate frequency of the linearized reduced system; while (b) sho ws tw o dominant frequencies close to the linear frequency resulting in a b eating phenomenon. (c) has E = 0 . 30, and is in the L ´ evy-w alk domain; the tra jectory transitions betw een differen t seemingly regular tra jectories with a p o wer-la w distribution. (d) has E = 1 . 87 and exhibits regular diffusion statistics. (e) E = 770 . 10 retains some regularity of motion, and (f ) E = 7 . 79 · 10 7 exhibits regular motion corresponding to the linearized high-energy system obtained by setting the rest length to zero. corresp onding to the ballistic motion c haracterizing the lo w energy , to the v alue 1, which characterizes the un- correlated random walks observed for the mo derate en- ergy v alues E s / 9 . E . 2 E s . The tra jectories in this regime displa y some regularity , follo wing a quasi-perio dic tra jectory with a constan t av eraged rotation velocity for a finite time, then transitioning to follo wing a different quasi-p eriodic tra jectory . This “sticking dynamics” [24], whic h results in the emergence of anomalous diffusion, is thus largely determined b y the regular quasi-p eriodic solutions of the Hamiltonian. The transitions themselves b et ween the seemingly quasi-p eriodic tra jectories, while unpredictable, are also related to the underlying phase space structure. Understanding the complex and subtle nature of this dynamics requires a deep understanding of the regular solutions of the system, and the correspond- ing structure of phase space. I II. LOW ENERGY REGIME In the low-energy regime E  E s , tra jectories exhibit a regular quasi-p eriodic motion. This regular motion seems to suggest that a p erturbativ e approach around the system’s zero-energy equilibrium would provide a go od description of the motion. How ever, as explained in 10 -2 10 0 10 2 10 4 10 6 Energy 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Largest Lyapunov exponent FIG. 3. Maximal Ly apunov exp onents as a function of en- ergy on a log scale, for v arious random initial conditions with zero angular momentum. At extremely low and high energies, the system b eha ves regularly and the maximal Ly apunov ex- p onen ts are zero. In mid-range energies, the system has a c haotic signature with a positive maximal Ly apunov exp onen t for most initial conditions. The Ly apunov exp onent calcula- tion w as done using the method describ ed in [25], using the Matlab program “Calculation Lyapuno v Exponents for ODE” v ersion 1.0.0.0 by V asiliy Gov orukhin. 4 [1], linearization in the Cartesian coordinates of Eq. (1) requires breaking the rotational symmetry of the prob- lem by c ho osing a sp ecific equilibrium p osition in the plane ab out which the linearization is performed. This t yp e of linearization, appearing in [26], conserves angu- lar momen tum only to leading order and fails to cap- ture the finite rotation of the triangle. Nev ertheless, the frequencies deriv ed from this linearization do match the strongest frequencies observed in the simulation, p 3 / 2 and √ 3. Capturing the true dynamics of the system in the low- energy regime requires a description of the system in its shap e subspace, as a deforming triangle instead of as three masses moving with a pairwise p oten tial. This pro cedure is p erformed in [1, 27–29] for different poten- tials and configurations of three-b ody systems, by a v ari- able change to three shape-space v ariables describing the shap e of the triangle and one angle determining the ori- en tation of the triangle in the plane. The c hoice of shape- space and orien tation co ordinates is a gauge c hoice that do es not impact the results. W e find that the most con- v enient choice for the shap e v ariables, presented in [28], is a Blo ch sphere represen tation of the tw o relative Ja- cobi co ordinates of the three masses, ρ 1 = p m 2 ( r 2 − r 1 ) and ρ 2 = q 2 m 3 ( r 3 − r 1 + r 2 2 ) (see Fig. 1). This v ari- able set is denoted w = ( w 1 , w 2 , w 3 ), and is giv en by w 1 = 1 2 ( ρ 2 1 − ρ 2 2 ), related to the isometric b end mo de; w 2 = ρ 1 · ρ 2 , related to the asymmetric b end mo de; and w 3 = ρ 1 ∧ ρ 2 , prop ortional to the oriented area of the triangle and related to the symmetric stretch mode. The rotation v ariable we use is θ , describing the angle be- t ween the line connecting m 1 and m 2 and the x axis. The cen ter-of-mass coordinates decouple from the rest of the system and are set to zero. The system’s inv ariance to rotations leads to the conserv ation of angular momentum J and allows one to deduce the orientational dynamics, expressed via ˙ θ , from the shap e-space dynamics through ˙ θ = J 2 w + w 2 ˙ w 3 − w 3 ˙ w 2 2 w ( w + w 1 ) , (2) where w = | w | . While ˙ θ is given by the explicit relation ab o ve as a function of w and ˙ w , one can sho w that θ can- not b e expressed as a function of w and ˙ w alone, which in turn enables the phenomenon of deformation induced rotation. The conserv ation of angular momentum th us yields a non-holonomic constraint for θ , and obtaining its v alue at a given time requires knowledge of the full dynamics of the system up to that time. As a result, θ b ecomes a sensitiv e measure for the system’s dynamics and correlations. P erforming a Routh reduction of the angular momen- tum [18], we set J = 0 to obtain the reduced shap e- space Hamiltonian describing the system’s zero angular momen tum motion, H red = w  p 2 1 + p 2 2 + p 2 3  + k 2 X ( r ij ( ~ w ) − L ) 2 , (3) where p i = ˙ w i 2 w , r ij ( ~ w ) = q 2  w − w · b ij  b 13 =  1 2 , √ 3 2 , 0  , b 12 = ( − 1 , 0 , 0) , b 23 =  1 2 , − √ 3 2 , 0  . The real space dynamics of the system can b e restored b y finding solutions to the reduced Hamiltonian (3) to find the shap e dynamics w ( t ) and substituting them into Equation (2) to obtain the orientation evolution θ ( t ). As pro ved in [18], ev ery solution of the original Hamiltonian (1) with J = 0 corresp onds to a solution of the reduced Hamiltonian (3), and vice v ersa, therefore it suffices to study solutions of (3). W e emphasize that the sim ulation results shown here are p erformed on the full, Cartesian system (1) while the p erturbativ e analysis is p erformed on the reduced system (3). A. P erturbation theory in the reduced shape space The reduction to shap e space allows us to expand the system not ab out a rest position but ab out its equilibrium shap e, the static equilateral triangle w 0 =  0 , 0 , mL 2 2  , p 0 = (0 , 0 , 0), thus allowing finite rotations of the triangle without breaking the small-p erturbation appro ximation. Redefining w i = w i 0 + α i ˜ w i , p i = α i − 1 ˜ p i for α 1 = α 2 =  2 L 4 m 3 3 k  1 / 4 , α 3 =  L 4 m 3 3 k  1 / 4 , and ex- panding in orders of  , we obtain the Hamiltonian as a p o wer series of the coordinates  ˜ w and  ˜ p . The first non- v anishing order is the linearized Hamiltonian, quadratic in the v ariables and th us integrable as a simple sum of harmonic oscillators. In action-angle v ariables it reads, H red =  2 2 E s r 3 2 ( I 1 + I 2 ) + √ 3 I 3 ! + O   3  , (4) where I j = 1 τ s E s ( ˜ w 2 j + ˜ p 2 j ) are the (dimensionless) action v ariables serving as generalized momenta, and their con- jugate co ordinates are the angle co ordinates denoted by φ i . The nonlinearity in the system is of geometric origin rather than constitutive, and in particular is not asso- ciated with an externally tunable expansion parameter; eac h of the individual springs is harmonic, and it is the geometric coupling of their strains that leads to non- linearit y . As a result the nonlinear effects increase con- comitan tly with the strains. The largest p ossible strain for a given total energy is b ounded and increases with the total energy . Thus, the total energy in the system can b e used to define an auxiliary expansion parameter,  ( E ), satisfying  (0) = 0 and monotonically increasing with the energy . Details of this rescaling are presented in the Supplementary Material (SM); in what follows we 5 use  as a dummy parameter in order to simplify nota- tion, recalling that rescaling can b e easily p erformed to yield a formal expansion for the p erturbative approac h in low energies. The linearized Hamiltonian (4) describ es three decou- pled harmonic oscillators corresp onding to the three vi- brational mo des of the planar triatomic molecule [30]: I 1 corresp onds to the asymmetric stretc h, I 2 to the b ending mo de and I 3 to the symmetric stretc h (see Fig. 1). W e note that the 1:1 resonance b etw een I 1 and I 2 is a re- sult of the symmetry of the system under consideration; c hanging, for example, one of the masses would remov e this frequency degeneracy . Substituting the solution of (4) in to the equation for ˙ θ 1 and a veraging out the fast os- cillations results in the following equation for the a verage angular velocity [28], ˙ θ 1 =  2 3 2 τ s p I 1 I 2 sin ( φ 2 − φ 1 ) . (5) As could b e inferred intuitiv ely , o verall rotation is a result of the phase difference φ 2 − φ 1 b et ween the tw o resonant oscillators I 1 and I 2 , the asymmetric stretch and the iso- metric b ending; the symmetric stretch osc illator I 3 has no rotational charge to first non-v anishing order around the equilibrium. In a comparison to simulations, we find that this ex- pression explains the o verall angular velocities well for lo w enough energies, but the fit deteriorates as the en- ergy is increased, see Fig. 4(c). The expansion to leading order also fails to account for the b eating phenomenon observ ed for some initial conditions (Fig. 2(b)). Seeking to impro ve the prediction for the rotation velocity as well as to account for the observ ed b eating one can attempt canonical p erturbation theory to higher orders; how ev er, due to the 1:1 resonance, the expansion diverges already at the next non-v anishing order. T o circumv ent this div ergence we recast the reduced Hamiltonian in Birkhoff normal form around its static equilibrium, using the metho d describ ed in [18, 31]. This impro ves the fit of the angular velocity and provides a go od description of the full dynamics observ ed in the reg- ular regime, see Fig. 5. The general pro cedure, presen ted in [31], is an iterative sc heme of wisely chosen canonical Lie transforms that puts the system in the form of a p olynomial series in action coordinates, where the series comm utes with its low est-order term. The main steps of the calculation of the normal form up to fourth or- der is presen ted in App endix B. In the main text we presen t the relev ant results, showing that the expansion to fourth order (presented in subsection B and C) pro- vides an accurate description of the basic phase space structure of the full system for low energies. F urther- more, the expansion to second order suffices to break the normal mo de frequency degeneracy which allows the application of Kolmogoro v-Arnold-Moser (KAM) theory (subsection D). B. 2 nd order Birkhoff normal form In order to recast the Hamiltonian in its Birkhoff nor- mal form to second order, we p erform a canonical change of co ordinates to a new set of v ariables, ( J, ψ ), which we will use throughout the rest of this section: J 1 = I 1 , J 2 = I 1 + I 2 , J 3 = I 3 and their conjugate angle co or- dinates ψ 1 = φ 1 − φ 2 , ψ 2 = φ 2 , ψ 3 = φ 3 . These v ari- ables naturally exhibit some of the in teresting b ehavior of the system, with J 3 as the energy contained in the area c hanges of the triangle, J 2 as the o verall energy contained in the resonant oscillators I 1 and I 2 , and J 1 as the en- ergy contained only in I 1 , satisfying J 1 ≤ J 2 . ψ 1 is the phase difference betw een the t w o resonant oscillators and describ es the av erage rotation direction of the triangle in the plane, with ψ 1 ∈ (0 , π ) and ψ 1 ∈ ( π , 2 π ) manifesting as counterclockwise and clo ckwise rotation, resp ectively . In these v ariables, the Birkhoff normal form of the sys- tem to second order is given by H (2) =  2 2 E s ( H 0 +  2 Z 2 ) (6) where H 0 = r 3 2 J 2 + √ 3 J 3 , Z 2 = 1 64  52 J 1 ( J 2 − J 1 ) sin 2 ψ 1 − J 2  5 J 2 + 6 √ 2 J 3  . (7) It is immediately apparent that the truncated system H (2) is integrable, with H 0 , J 2 and J 3 conserv ed quan- tities. Therefore, in order to visualize the dynamics of the truncated system, it suffices to consider the tw o- dimensional phase plane J 1 J 2 - ψ 1 giv en v alues J 3 ≥ 0 and J 2 > 0 1 , as seen in Fig. 4. This phase plane has the structure of a finite cylinder, with 0 ≤ J 1 /J 2 ≤ 1 and ψ 1 ∈ [0 , 2 π ] an angle v ariable, and can b e un- folded onto the plane (see Fig. 4). The system has n ullclines at J 1 /J 2 = 0 and J 1 /J 2 = 1, fixed lines at ψ 1 = k π for k ∈ N and t wo distinct elliptic fixed points at ( J 1 /J 2 , ψ 1 ) = (1 / 2 , π / 2) and ( J 1 /J 2 , ψ 1 ) = (1 / 2 , 3 π / 2). Solutions of H (2) p erform closed orbits around the dis- tinct fixed p oints and cannot cross the rectangles drawn b y the nullclines and the fixed lines. These tra jecto- ries display a b eating phenomenon where energy p eri- o dically transfers b et ween the resonant oscillators, with more energy contained in I 1 ( I 2 ) when J 1 /J 2 < 1 / 2 ( J 1 /J 2 > 1 / 2). The integrable dynamics of the truncated system lie in the shap e space of the system, and can b e pulled back to obtain the corresp onding rotation of the triangle in 1 When J 2 = 0, J 1 v anishes as well, corresp onding to equilateral triangles. This is an in tegrable family of sp ecial symmetries, with no discernible impact on the observable phase space when J 2 > 0. 6 0 50 100 150 200 t[ s ] -15 -10 -5 0 5 10 15 (t) [rad] 0 0.05 0.1 0.15 0.2 0.25 0.3 Energy [arbitrary units] 0 0.005 0.01 0.015 0.02 [rad/ s ] Full system Linear approximation Birkhoff approximation D ˙ ✓ E 0 50 100 150 200 t[ s ] -15 -10 -5 0 5 10 15 (t) [rad] 0 0.05 0.1 0.15 0.2 Energy [arbitrary units] 0 1 2 3 4 5 6 7 10 -3 Full system Linear approximation Birkhoff approximation (a) (b) (c) (d) FIG. 4. (a), (b) and (c) sho w the dynamics of the truncated system H (2) . (a) shows typical tra jectories of H (2) pro jected onto the ψ 1 vs. J 1 /J 2 plane. (b) displa ys the corresp onding rotations of these tra jectories. (c) is a long-exp osure image of their dynamics of three out of the four tra jectories, with each mass colored in a different color. (d) sho ws a comparison of the full system’s av erage rotation velocity as a function of the energy with the linear approximation and the muc h improv ed Birkhoff 2nd-order approximation given b y H (2) . the plane. As seen in Fig. 4, rotation around the fixed p oin t ψ 1 = π / 2 manifests as a negativ e angular v elo cit y , while rotation around the other fixed p oint ψ 1 = 3 π / 2 as p ositive angular v elo city . T ra jectories that are close to the fixed p oints hav e a smaller b eating frequency and a larger absolute angular velocity than tra jectories that pass closer to the rectangles’ b orders, in which the b eat- ing is v ery apparent, for the same energy . Despite the b eating, since the motion around the fixed p oints is p e- rio dic, the av eraged rotational velocity of a given tra jec- tory is constant in all cases. C. Birkhoff normal form to higher orders Of course, a priori there is no guarantee that the nor- mal form should describ e or ev en approximate the dy- namics of the full system that we observe in simulations. The normal form series do es not necessarily con verge, and the theory is guaran teed to hold only for some neigh- b orhoo d of the static equilibrium configuration. W e can c heck compatibility by comparing the simulated dynam- ics of the full system (1) with the predicted normal-form dynamics of H (2) , pro jecting the tra jectories on to the unfolded J 1 /J 2 - ψ 1 plane. As can be seen in Fig. 5b, the tw o elliptic fixed p oin ts of H (2) are clearly visible in the full system. T ra jectories of the full system that initialize close enough to the fixed p oints shadow the truncated system’s tra jectories, with similar b eating fre- quencies and a similar ov erall angular velocity of the ro- tating triangle (Fig. 4c). Ho wev er, despite this high compatibilit y , some prominent elements of the dynam- ics are not captured by H (2) . The lines ψ 1 = k π for in teger k , are not fixed for the full system; rather, sim u- lations of the full syste m sho w an elliptic fixed p oin t at ( J 1 , ψ 1 ) = ( J 2 / 4 , π ), around which b eating tra jectories with zero ov erall angular velocity are distinctly apparent in the simulations. In order to capture these features, w e now come to con- sider the Birkhoff normal form to the next non-v anishing order, given by H (4) =  2 2 E s ( H 0 +  2 Z 2 +  4 Z 4 ) , (8) where Z 4 = a 0 + J 1 ( a 1 + b 1 cos 2 ψ 1 ) + J 2 1 ( a 2 + b 2 cos 2 ψ 1 ) + J 3 1 ( a 3 + b 3 cos 2 ψ 1 ) . (9) The co efficients { a i , b i } 3 i =0 are functions of J 2 and J 3 , presen ted in full in App endix 3 along with the full cal- culation. H (4) still conserves J 2 and J 3 , retaining the in tegrability of the normal form. The distinct elliptic fixed p oints of H (2) are also fixed p oin ts in the phase space of H (4) . How ev er, the ψ = k π , k ∈ Z lines lose their stability; instead, tw o new fixed p oin ts e merge on eac h line, an elliptic fixed p oin t at ( J 1 , ψ 1 ) =  J 2 4 , k π  and a h yp erb olic fixed p oin t at ( J 1 , ψ 1 ) =  3 J 2 4 , k π  . Also, four h yp erb olic fixed p oin ts app ear on the nullcline of J 1 = 0; see Appendix C for the full fixed p oint analysis. The fixed p oints separate the tra jectories into tw o classes of tra jectories, those en- circling the elliptic fixed p oin ts and those migrating along the nullclines of J 1 = 0 and J 1 = J 2 . 7 (a) (b) 0 0.25 0.5 0.75 1 J 1 /J 2 0 /2 3 /2 2 FIG. 5. (a) Phase space of H (4) compared with (b) Poincar ´ e sections of the full system (1) at a low energy E = 0 . 0075. (a) is divided into four dynamical regions: the purple region, in which tra jectories encircle the fixed p oin t J 1 = J 2 / 2 , ψ 1 = 3 π / 2, describing clo c kwise rotation of the triangle; the blue region of tra jectories encircling J 1 = J 2 / 2 , ψ 1 = π/ 2 describing coun terclo ckwise rotation; the orange region of tra jectories encircling J 1 = J 2 / 4 , ψ 1 = 0; and the green region of tra jectories migrating along the nullcline J 1 = J 2 . (b) shows that for low enough energies, the full system follows the regular structure of H (4) to a high level of accuracy , with six t ypical tra jectories shown in the purple region, five in the blue region, three in the orange region and one in the green region.. Pulling bac k from shap e space to the space of triangles in the plane, the rotation resulting from the dynamics of H (4) are similar to that of H (2) close enough to the fixed p oin ts shared by the systems. Ho w ever, the dynamics are different around the new phase space features: tra- jectories going around ( J 2 / 4 , k π ) p erio dically rotate in b oth directions in real space, with an o verall v anishing a veraged rotation rate. The tra jectories follo wing the n ullclines also do not p erform ov erall rotation. The het- ro clinic tra jectories b et ween the h yp erb olic fixed p oints define the b oundary b et ween the different dynamical re- gions. As can b e seen in Fig. 5, these dynamics are indeed compatible with the full system to a high degree. P oincar´ e sections of the full dynamics pro jected onto the ( J 1 /J 2 , ψ 1 ) plane reveal exactly the same fixed p oints as calculated from H (4) . D. Lifting the frequency degeneracy The harmonic three mass system sho ws a strong p er- sistence of regular solutions for a large range of ener- gies. This suggests the applicability of the KAM theo- rem, which guaran tees p ersistence of most quasi-p erio dic orbits in almost-integrable systems if the integrable part satisfies some non-degeneracy frequency condition. Un- fortunately , the harmonic expansion (to low est order), H 0 , shows a one-to-one resonance, and thus cannot serve as the base of a KAM expansion. Therefore, to express our system as an almost-in tegrable system, we write the full Hamiltonian as H = H (2) + ( H − H (2) ). The Birkhoff expansion, detailed in the SM, shows that the remain- der in the parentheses is a p o wer series, H − H (2) = P ∞ n =2 P ( n ) ( J , ψ ) where P ( n ) is a monomial of order n in the action co ordinates. The auxiliary expansion pa- rameter  ( E ), a monotonically decreasing function of the maximal p ossible stretc h given the energy E (see SM), can then b e used to rescale the action co ordinates, so that 0 ≤ J i ≤ 1 for i = 1 , 2 , 3. Thus, the full Hamilto- nian is rewritten as a p ow er series in  ( E ) multiplying terms of order 1. Hence, so long as the en tire remainder is small, w e can treat our system as an almost-integrable system, considering ( H − H (2) ) as the perturbation to the in tegrable and non-degenerate H (2) . The KAM theorem do es not provide a realistic b ound on what consists a small enough p erturbation for the- orem to b e applicable. Nev ertheless, we can c heck the frequency conditions required for the theorem by p er- forming a canonical change of v ariables to action-angle 8 v ariables and calculating the frequencies:x H (2) =  2 2 E s r 3 2 J 2 + √ 3 J 3 −  2 64  5 J 2 2 − 6 √ 2 J 2 J 3 + 13 K 2 1   ; (10) ω 1 =  2 13 32 K 1 , ω 2 = r 3 2 −  2 32  5 − 3 √ 2  J 2 + 13 K 1  , ω 3 = √ 3 +  2 3 √ 2 32 J 2 ; (11) where K 1 = 2 p J 1 ( J 2 − J 1 ) sin ψ 1 is the new conserved quan tity emerging from the in tegrable system. Note that it is prop ortional to the linear slop e prediction Eq. (5); indeed, the sign of K 1 indicates the ov erall direction of rotation of the triangle, see Fig. 6(c). The new frequencies asso ciated with the angle co ordi- nates ha ve corrections of order  2 whic h dep end on the action co ordinates, th us remo ving the degeneracy of the linearized system. It is easy to chec k that both the non- degeneracy and the iso energetic non-degeneracy condi- tions stated in the KAM theorem are satisfied for small enough v alues of J 2 , J 3 . Under these conditions, the KAM theorem assures that most in tegrable tori p ersist under small p erturbations to the Hamiltonian for an y en- ergy v alue that is small enough. T aking into account H (4) as the integrable part would add corrections of order  4 , retaining this degeneracy lifting. The loss of integrabilit y is exp ected to manifest first around resonant tori, ov ertaking most of the phase space gradually as the perturbation gro ws. This picture is com- patible with our n umerical exp eriments and provides a p ossible explanation for the go o d fit b et ween the trun- cated and the full system’s dynamics for lo w enough en- ergies. F or short times, this shadowing of the tra jecto- ries of the truncated in tegrable Hamiltonian by the full Hamiltonian tra jectories remains as the energy is further increased, as we show next. E. The L´ evy W alk regime A t energies in the range E s / 15 . E . E s / 9 most tra- jectories are no longer regular: chaotic dynamics charac- terized by a p ositive Lyapuno v exp onen t inhibit most of phase space. A t the early stages of this regime, the cor- resp onding rotational dynamics resemble the sto c hastic L ´ evy-w alk model [1], with bouts of constant angular ve- lo cit y interrupted by abrupt orien tation reversal even ts. In Fig. 6 we plot the angular dynamics in this regime, alongside a pro jection of phase space onto the ( J 1 /J 2 , ψ 1 ) plane, where the system is shown to follo w the integrable structure describ ed b y H (4) . The pro jection indicates that the observed tra jectories migrate betw een the differ- en t fixed p oin ts of H (4) , sticking to oscillatory tra jectories around each of the stable fixed p oints for long times. W e further observ e that the transitions b etw een the distinct neigh b orho o ds of the fixed p oin ts o ccurs near the saddle p oin ts lo cated at J 1 /J 2 = 3 / 4 , ψ = π k for k ∈ Z , and the transition times ob ey a p ow er law distribution. A t the lo west energies in which we observe L´ evy- w alks each b out b et ween orientation reversal even ts b ears great resemblance to the corresp onding regular tra jec- tory around the same fixed p oint, to the extent that it is difficult to differentiate b et ween regular and L ´ evy-w alk tra jectories just by examining them for short times in b et ween transitions. As the energy is increased the tran- sitions b et ween the neigh b orho ods of the fixed p oin ts b e- come more frequent and occur ov er an increasingly wider region. As the energy approaches E s / 9 from b elo w it seems that there is no longer an y barrier separating the basins of the distinct fixed p oints, and the b outs grad- ually lose their coherence and similarity to the regular solutions. Nonetheless, the squared angular mean dis- placemen t still ob eys fractional statistics [1]. In lo w-dimensional systems, d ≤ 2, the emergence of p ow er-law statistics is well-understoo d, attributed to the breakdo wn of Kolmogoro v-Arnold-Moser (KAM) tori creating partial transp ort barriers which can b e crossed b y chaotic tra jectories at a slow rate in a phenomenon commonly referred to as sticking or trapping [5 – 10]. Ho wev er, despite the robustness of this phenomenon [4, 11 – 14], a general framework for the origin of p ow er- la w statistics in high-dimensional mixed Hamiltonian sys- tems con tin ues to elude current understanding [5, 15–17]. In the three-b o dy harmonic system, the coherent b outs creating the p ow er-law statistics strongly resem ble their lo wer-energy regular counterparts, indicating a partial trapping of chaotic tra jectories around regular islands for finite times. A quantitativ e analysis of the dynamical mec hanism b ehind this phenomenon would require a deep understanding of the regular b ehavior of the non-linear system. In mixed Hamiltonian systems, fractional statistics are ubiquitous [4, 11 – 14]. In systems with t wo degrees of free- dom, where regular tori create barriers in phase space, the origin of these anomalous statistics is well under- sto od. Generally , as the KAM tori break up, they lea ve in their wak e a hierarchical structure of smaller tori that create partial barriers of transp ort. Chaotic tra jectories can cross these barriers, but this typically takes a long time, resulting in the fractional statistics [7, 17]. Ho wev er in higher dimensional systems, the phase space mec h- anism creating and controlling the observed fractional statistics is not yet fully understo o d [5]. Although the KAM tori break up in a similar manner, they no longer separate phase space into imp enetrable regimes. Th us c haotic tra jectories can theoretically get as close as they lik e to the surviving tori. In these systems, for any p er- turbation strength, the phase space is connected b y a web of resonant c hannels known as the Arnold web surround- 9 ing the sufficiently non-resonant KAM tori. Action v ari- ables can drift along these channels in a pro cess known as Arnold diffusion and thus transition from the neigh- b orhoo d of one surviving torus to another. In our sys- tem, the great resem blance of the low-energy L´ evy-w alk tra jectories to regular solutions and the narro w channel of transfer are reminiscent of the Arnold diffusion phe- nomenon. On the other hand, as the energy rises and the transition region gro ws, the H 4 phase space structure loses its coherence and the p o wer-la w statistics seem to originate from a partial trapping of tra jectories around the regular fixed p oin ts. A combination of the t wo phe- nomenon could explain the surprising phenomenon of a gradual, seemingly contin uous decrease of the anomalous exp onen t from the ballistic to the random w alk regime as the energy grows, as observed in [1], as opp osed to the single anomalous exp onen t found in [10, 32]. This w ork provides the backbone that would b e required for a systematic study of these concepts, by iden tifying the un- derlying almost-integrable appro ximation controlling the dynamics in the transition of the full system from regu- lar b ehavior to c haos. These allow a calculation of the KAM tori and the surrounding Arnold web. A quantita- tiv e study of these ideas is left to future w ork. IV. HIGH ENER GY REGIME As the energy is increased b ey ond E s / 15, the regu- lar structure gradually disapp ears. In the range E s / 9 . E . 2 E s , almost all tra jectories are observed to cov er the en tirety of phase space, and the statistics resemble an uncorrelated random walk [1]. Ho wev er, at E ≈ 2 E s a single frequency b egins to dom- inate the dynamics and the system app ears to approac h regularit y again. This apparent regularit y may b e easily explained by observing that for extremely high energies, the rest length is effectively forgotten. Th us we ma y exp ect the system to resem ble the in tegrable harmonic three-mass system with zero rest lengths [33] at high enough energies. In this regime the reduced system (3) is less instructive, since its normal modes do not coincide with those of the zero rest length harmonic three-mass problem. Therefore, we compare the observed dynamics with the dynamics of the full Cartesian Hamiltonian (1) with zero rest lengths: H = 3 X i =1 p 2 i 2 m + X k 2 r 2 ij , (12) whic h displa ys the twice-degenerate linear frequency √ 3 in units of 1 /τ s . As sho wn in Fig. 7o, the solutions of (12) are in ex- cellen t agreement with the simulation results for high enough energies. Similarly to the low energy regime, the regular solution of the high energy regime displays a constant a verage angular slop e. How ever, unlike the solution for low energies, this slop e is comprised of steps: 0 0.25 0.5 0.75 1 J 1 /J 2 0 /2 3 /2 2 1 (a) (b) 0 1000 2000 3000 t[ s ] -10 -5 0 (t) [rad] 0 1000 2000 3000 t[ s ] -0.3 -0.2 -0.1 0 0.1 0.2 K 1 [energy] (c) FIG. 6. (a) A typical orientational tra jectory at energy E = 0 . 381. (b) Pro jection of the phase space dynamics onto a Poincar ´ e section of the ( J 1 /J 2 , ψ 1 ) plane. The p ow er-law statistics observ ed in this energy regime correspond to a stic k- ing of the irregular tra jectories close to the approximated in- tegrable system’s fixed p oints for long times. (c) The action v ariable K 1 emerging from H 2 is shown in light blue, with its mo ving av erage h K 1 i on top, av eraged ov er 100 time units. h K 1 i < − 0 . 03 is colored in blue and corresponds to a descend- ing rotation angle; h K 1 i > 0 . 03 is in magen ta and corresp onds to an increasing rotation angle; and − 0 . 03 < h K 1 i < 0 . 03 is in yello w and corresp onds to zero av eraged rotation. K 1 do es not differentiate b etw een tra jectories that encircle the J 1 /J 2 = 0 . 25 (orange in (a,b)) and tra jectories that trav el along the J 1 /J 2 = 1 b order (green in (a,b)), but may identify a transition b et ween them. It is obvious from (c) that K 1 is not a conserved quantit y of the full Hamiltonian, nor is it a constant along seemingly ballistic b outs. Nonetheless, its mo ving av erage is indicative of the different regimes. discrete angular increment ev ents in b et ween which the angle remains approximately constant. This feature may b e explained by observing that for high energies, at ev ery oscillation the three-mass triangle undergo es tw o orien- tation reversal even ts. Both orientation rev ersal even ts o ccurring in a single oscillation increase (or decrease, de- p ending on initial conditions) θ ( t ) by π . Hence the ex- p ected a veraged slop e has a constant v alue of √ 3 ≈ 1 . 73, as observed. The step structure is preserved at mo derately high en- ergies, where the rest lengths are not completely negligi- ble, and is observed even at energies that are v ery close to the energy scale, E & 2 E s , see Fig. 7 (l,m,n). F or these mo derate energies, the dominant frequency drifts a wa y from √ 3 and the sp ectrum fills up (Fig. 7 b,c,d), and for low enough energies θ ( t ) changes its av erage rotation direction in a manner resembling the L´ evy walk region. A statistical and p erturbativ e analysis of the Hamilto- nian in this energy regime is left to future studies, but it is clear that a similar approach to the p erturbative tec hniques used in the lo w energy regime may b e useful 10 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) p 3 t [ ⌧ s ] ✓ ( t ) P ( ! ) ! [1 / ⌧ s ] linear theory full system FIG. 7. (a-e) Po wer sp ectrum of the x co ordinate of m 1 for v arious typical tra jectories with rising energies v alues from E ∼ 1 to E ∼ 10 6 . (f-j) Long-exp osure images of the corre- sp onding tra jectories. (k-o) Orientation as a function of time of the full system (blue) and the linear prediction obtained b y solving the zero-rest length system with the same initial conditions (green). As the energy gets larger, the high-energy limit is a b etter approximation to the full system dynamics. The energy v alues of the different plots are (a,f,k) E = 1 . 47, (b,g,l) E = 9 . 51, (c,h,m) E = 95 . 42, (d,i,n) E = 1354, (e,j,o) E = 2 . 11 × 10 6 . in analyzing the approach of the system to high energies where the system b ehav es like its linear appro ximation. V. SUMMAR Y AND DISCUSSION Despite the apparent simplicit y of the harmonic three- mass system with finite rest lengths, the system displa ys a rich v ariety of dynamics, controlled mainly by the sys- tem’s ov erall energy . F or very low energies the system displa ys constant angular velocity rotation with zero an- gular momen tum, while for mo derate energies chaos en- sues, indicated by a p ositive Lyapuno v exp onent, and an orientational random walk is measured. Gradually increasing the energy from very low to mo derate v alues rev eals statistics of an orientational L´ evy walk, in which the exp onent α contin uously v aries with the energy , in- terp olating b et ween the v alues α = 2 (ballistic rotation) and α = 1 (random walk). F urther increasing the energy b ey ond the random walk region, the system gradually “forgets” its finite rest lengths and the systems tra jecto- ries regain regularity . In the chaotic regimes, while the observed tra jectories share many characteristics with regular solutions of in- tegrable approximations of the Hamiltonian, no explicit solutions are av ailable. Th us identifying the fixed p oin ts of the system’s Hamiltonian, the structure of the regu- lar solutions and the geometry of the phase space for the in tegrable approximations of the Hamiltonian is key to understanding the exotic phenomena that the full system displa ys. F urthermore, observing the ric h phenomena that the harmonic three m ass system displa ys requires to sim ulate the underlying chaotic Hamiltonian system to very long times ( ∼ 10 7 Ly apunov times). Such a task is not commonly carried out, primarily b ecause of the difficult y in interpreting the result and identifying the real system it describ es [34]. Understanding analytically the building blo c ks from which the differen t parts of the full tra jectory is comp osed not only allo ws to understand the origin of the observed L ´ evy walks, but also serves to cast meaning to the observed tra jectories as typical mem b ers in a collection of statistically similar tra jecto- ries that cov er the chaotic comp onen t of phase space. T o these ends, in this work we mapp ed and characterized the fixed p oints and regular solutions of b oth the lo w and high energy integrable appro ximations of the Hamil- tonian of the harmonic three mass system. V ery low energy tra jectories displa y a constant angu- lar velocity rotation with zero angular momentum that is w ell captured b y reducing the system to its intrinsic (shap e) space and linearizing the system ab out its equi- librium shap e. As the energy is increased non-linear ef- fects including lifting of the frequency degeneracy and b eating dramatically change the v alue of the constan t angular velocity . Capturing these v ariations required a p erturbativ e approach. Due to the 1:1 resonance in the system’s linearization, canonical p erturbation theory di- v erges, therefore w e utilized the Birkhoff normal form ex- pansion. Expanding the Hamiltonian in Birkhoff normal forms to 4th order yielded a v ery go o d agreement with the observed av erage angular velocity , and was shown to accurately predict the phase space structure observed for the full system. The action angle v ariables inherited from these approximations capture the b ehavior of the system near the onset of chaos. F or higher energy , when L ´ evy walks b ecome more pronounced these action v ari- ables are no longer constan t even along seemingly ballis- tic bouts. Nonetheless, the predicted phase space struc- ture still underlies the full dynamics, and the L´ evy walks can b e decomp osed to b outs that dw ell in the vicinity of regular tra jectories rotating clo c kwise and anticlockwise at a constant pace, and p o wer-la w distributed transition b et ween these tra jectories. As the system approaches regular random w alk with α = 1 the structure of phase space predicted from the in tegrable approximations ceases to describ e the system. F or a narrow strip of energies in the random w alk re- gion the system app ears lacking an underlying structure. Ho wev er, as the energy of the system is further increased 11 a new structure emerges. When the t ypical mass sepa- ration significantly exceeds the rest length, L  h| r ij |i , the rest lengths are effectively lost, and for 10 4 E s . E the observ ed tra jectories seem regular again, and can b e explained considering the system with v anishing rest lengths. This new structure b egins to b e apparent al- ready at energies E s < E as a single frequency starts to dominate the p o wer sp ectrum of the system, while still in the chaotic regime. F or energies in the intermediate range E s / 15 . E . E s / 9, anomalous p ow er-law statistics of the system’s v ariables are measured [1]. The phase-space mec hanism b ehind p o wer-la w correlations and corresponding anoma- lous diffusion of measurables in systems with a high phase space dimension is not w ell-understo o d, and ma y be at- tributed to Arnold diffusion, stickiness or some combina- tion of the tw o [5]. The quantitativ e understanding of the regular structure achiev ed in this work is crucial in order to understand and quantify the anomalous region, and to differen tiate b etw een the mechanisms responsible. In our system we observe that at the onset of the region, tra- jectories sp end long times circling one of the lo w-energy Birkhoff expansion fixed p oin ts, resembling their corre- sp onding regular tra jectories, b efore transitioning to a differen t fixed p oint through a narrow transfer channel around a hyperb olic fixed p oin t. This scenario is rem- iniscen t of an Arnold diffusion mechanism. As the en- ergy rises, the KAM islands shrink and this description gradually loses its coherence, resulting in an anomalous exp onen t that app ears to in terp olate smo othly b et ween ballistic and random-walk v alues [1]. While we presently cannot pro ve so, we b elieve b oth Arnold diffusion and stic ky dynamics dominate the system’s b ehavior for dif- feren t energies, partially explaining the smooth interpo- lation b et ween the ballistic and regular diffusiv e regimes. ************* APPENDICES A. Obtaining a small parameter  ( E ) In formulating the system using p erturbation theory ,  was an auxiliary v ariable that only served as a dummy parameter to ease the expansion. Ho wev er, using the ge- ometrical constraints of the system we can provide an estimate for ε ( E ), and use it to rescale the parameters so that they remain b ounded. This is a sp ecial feature of the spring-mass system, as opp osed to some other c haotic systems suc h as the three-b o dy gravitational problem: since the full Hamiltonian is p ositive-definite in its pa- rameters, the o verall energy of the system limits the ki- netic energy that the masses can gain, and b ecause of the geometry the masses cannot drift farther a wa y from eac h other than a certain radius. This allo ws us to place b ounds on the action v ariables given an energy v alue E , 0 ≤ I 1 ≤ 4 q 2 3 E 3 k  L − √ 2 q E k  2 , 0 ≤ I 2 ≤ 4 q 2 3 E 3 k  L − √ 2 q E k  2 , 0 ≤ I 3 ≤ 4 E 3 √ 3 k  L − q 2 E k  2 W e thus define  ( E ) to be the larger of the three: ε ( E ) 2 = 4 q 2 3 E 3 k  L − √ 2 q E k  2 , (13) defined so that at a giv en energy E the action v ariables I j cannot surpass ε ( E ) 2 . An estimate of the energy at whic h p erturbation theory is exp ected to break down is giv en by comparing ε ( E ) to 1, o ccurring at E ≈ 0 . 66. Indeed as n umeric show, this v alue is close to the energy at which we see an onset of c haos. F urther, by rescaling the action parameters I j = ε ( E ) 2 ˜ I j , we know that their range is alwa ys 0 ≤ ˜ I j ≤ 1, and ε ( E ) is a monotonically increasing function of E , satisfying  ( E = 0) = 0. Therefore, for small enough en- ergies the bulk of the energy is contained in low orders of the ε expansion, constraining the remainder and provid- ing further justification of the applicabilit y of p erturba- tion theory techniques to analyze the system as nearly- in tegrable. B. Birkhoff Normal F orm to 6th Order Giv en an m-dimensional Hamiltonian H with an el- liptic fixed p oint at the origin, consider the linearized 12 Hamiltonian ab out its fixed p oin t, H 0 = P m i =1 ω i p 2 i + q 2 i 2 . Then the Birkhoff normal form theorem states that for an y p ositiv e integer N ≥ 0 there exists a neighborho o d U N ⊂ R 2 n of the origin and a canonical transformation T N : U N → R 2 n that brings the full system to its Birkhoff normal form up to order N: H ( N ) := H ◦ T N = H 0 + Z ( N ) + R ( N ) (14) where Z ( N ) is a p olynomial of degree N+2 that Pois- son comm utes with its leading order expansion ab out the fixed p oin t, H 0 , i.e.  H 0 , Z ( N )  ≡ 0 , and R ( N ) is small, i.e.   R ( N ) ( x )   ≤ C N | x | N +3 , ∀ x ∈ U N . A pro of of this theorem is given in [31]. It is a con- structiv e pro of with a general recip e for obtaining the Birkhoff normal form up to any desired order N ∈ N , giv en a Hamiltonian with an elliptic fixed p oint at the origin. Here we presen t the main steps of the construc- tion for our system. The recip e is based on a series of Lie co ordinate trans- forms chosen such that the p olynomial correction Z ( N ) P oisson commutes with H 0 . A Lie transform of co or- dinates is a canonical c hange of v ariables induced by some generating function χ . Assume we hav e a p oly- nomial g ( p, q ) of order n + 2, and a Lie transform gener- ator χ ( p, q ), which is a p olynomial of order m . Consider φ t χ = ( p ( t ) , q ( t )), the prop ogation of the v ariables p and q according to a Hamiltonian given by χ . W e seek to express the original p olynomial g ( p, q ) estimated at the prop ogated co ordinates: g ( p ( t ) , q ( t )) ≡ g ◦ φ t χ . Setting t = 1, the new p olynomial can b e written as a p ow er series in the order of the polynomials, g ◦ φ 1 χ = X k ≥ 0 g k (15) where g 0 := g , g k = 1 k { χ, g k − 1 } , k ≥ 1 (16) and the order of the p olynomial g k is n + k m . Consider now a p olynomial Hamiltonian expanded in p o wers of the co ordinates and momenta about its elliptic fixed point, H =  2 H 0 ( p, q ) + P ∞ n =1  n +2 P n ( p, q ), where P n ( p, q ) is a sum of monomials of order n + 2, of the form q L p n +2 − L . This Hamiltonian is already in Birkhoff nor- mal form to zeroth order. F or any first-order p olynomial χ 1 , the corresp onding Lie transform of H leads to the ordered form: H ◦ φ χ 1 =  2 H 0 +  3 ( P 1 + { χ 1 , H 0 } )+  4 ( P 2 + { χ 1 , P 1 } + { χ 1 , { χ 1 , H 0 }} ) + O (  5 ) (17) As P 1 + { χ 1 , H 0 } is a p olynomial of order 3, choosing χ 1 suc h that this term commutes with H 0 will bring Eq. (17) to its Birkhoff normal form up to 1st order. In general, obtaining an n ’th degree Birkhoff normal form is done iterativ ely . Consider a Hamiltonian giv en in its Birkhoff normal form up to order n − 1, i.e. H ◦ T n − 1 =  2 H 0 + Z ( n − 1) + R ( n − 1) : Z ( n − 1) is a p olynomial of order n − 1 that commutes with H 0 , and R ( n − 1) is of order ≥ n . W riting the remainder R ( n − 1) as a series of monomials of increasing order, R ( n − 1) = P ∞ k = n  k R k , a Lie transform induced b y a generating p olynomial χ n of order n + 2 will result in the following form for the Hamiltonian: ( H◦T n − 1 ) ◦ φ χ n =  2 H 0 + Z ( n − 1) +  n ( { χ n , H 0 } + R n )+ O ( n +1) . (18) Then, χ n is c hosen suc h that { χ n , H 0 } + R n P oisson com- m utes with H 0 ; the remaining terms will b e of higher orders from the construction. In particular, this implies that the Birkhoff normal form to 2nd order is obtained by c ho osing the 4th de- gree p olynomial χ 2 suc h that { χ 2 , H 0 } + R 2 comm utes with H 0 , where R 2 ≡ P 2 + { χ 1 , P 1 } + { χ 1 , { χ 1 , H 0 }} . F or further details, including the metho d used to choose the functions χ k , see [31], which includes a result about the time-scales at whic h the truncated system H 0 + Z ( n ) ma y b e considered instead of the full system. F ollowing this recip e, we obtain the follo wing trun- cated Birkhoff normal form of our system to order 4: H (4) =  2 H 0 +  4 Z 2 +  6 Z 4 H 0 = r 3 2 J 2 + √ 3 J 3 , Z 2 = − 1 64  52 J 1 ( J 1 − J 2 ) sin 2 ψ 1 + J 2  5 J 2 + 6 √ 2 J 3  Z 4 = a 0 + J 1 ( a 1 + b 1 cos 2 ψ 1 ) + J 2 1 ( a 2 + b 2 cos 2 ψ 1 ) + J 3 1 ( a 3 + b 3 cos 2 ψ 1 ) , (19) where J 1 = I 1 , J 2 = I 1 + I 2 , J 3 = I 3 , ψ 1 = φ 1 − φ 2 , ψ 2 = φ 2 , ψ 3 = φ 3 , and { I k , φ k } are the action-angle v ariables asso ciated with the linearized Hamiltonian H 0 , I k = 1 τ s E s ( ˜ w 2 k + ˜ p 2 k ), and: a 0 = 4606 √ 2 J 3 2 − 12401 J 2 2 J 3 − 38752 √ 2 J 2 J 2 3 + 8736 J 3 3 344064 √ 3 , a 1 = − J 2  156017 √ 2 J 2 + 1674 J 3  344064 √ 3 , b 1 = J 2  837 J 3 − 43505 √ 2 J 2  172032 √ 3 , a 2 = 199522 √ 6 J 2 + 837 √ 3 J 3 516096 , b 2 = 124514 √ 6 J 2 − 837 √ 3 J 3 516096 , a 3 = − 45005 28672 √ 6 , b 3 = − 9001 q 3 2 28672 . (20) 13 [1] O. S. Katz and E. Efrati, Physical review letters 122 , 024102 (2019). [2] T. Geisel, A. Zacherl, and G. Radons, Ph ysical review letters 59 , 2503 (1987). [3] T. Geisel, in L´ evy flights and r elate d topics in physics (Springer, 1992) pp. 151–173. [4] V. Zaburdaev, S. Denisov, and J. Klafter, Reviews of Mo dern Ph ysics 87 , 483 (2015). [5] S. Lange, A. B¨ ack er, and R. Ketzmeric k, EPL (Euro- ph ysics Letters) 116 , 30002 (2016). [6] R. S. MacKa y , J. D. Meiss, and I. C. P erciv al, Ph ysica D: Nonlinear Phenomena 13 , 55 (1984). [7] J. D. Meiss and E. Ott, Physica D: Nonlinear Phenomena 20 , 387 (1986). [8] G. Cristadoro and R. Ketzmeric k, Physical review letters 100 , 184101 (2008). [9] O. Alus, S. Fishman, and J. D. Meiss, Physical Review E 90 , 062923 (2014). [10] O. Alus, S. Fishman, and J. D. Meiss, Physical Review E 96 , 032204 (2017). [11] V. Latora, A. Rapisarda, and S. Ruffo, Ph ysical Review Letters 83 , 2104 (1999). [12] F. Cagnetta, G. Gonnella, A. Mossa, and S. Ruffo, EPL (Europh ysics Letters) 111 , 10002 (2015). [13] M. Shlesinger, B. W est, and J. Klafter, Ph ysical Review Letters 58 , 1100 (1987). [14] G. M. Zaslavsky , Physics rep orts 371 , 461 (2002). [15] C. Danieli, D. Campb ell, and S. Flach, Physical Review E 95 , 060202 (2017). [16] S. Das and A. B¨ ac ker, arXiv preprin t (2019). [17] J. D. Meiss, Chaos: An Interdisciplinary Journal of Non- linear Science 25 , 097602 (2015). [18] V. I. Arnold, V. Kozlov, and A. Neishtadt, “Dynamical systems iii. encyclop edia of mathematics, vol. 3,” (1988). [19] A. Guichardet, Ann. Inst. H. Poincar ´ e 40 , 329 (1984). [20] F. Wilczek and A. Shap ere, Ge ometric phases in physics , V ol. 5 (W orld Scientific, 1989). [21] R. Montgomery , Fields Inst. Commun 1 , 193 (1993). [22] F. J. Beron-V era, Symplectic in tegrators (h ttps://www.mathw orks.com/matlab central/fileexc hange/7686- symplectic-in tegrators), MA TLAB Central File Ex- c hange. Retrieved 2016. [23] E. Hairer and G. W anner, in Encyclop e dia of Applie d and Computational Mathematics (Springer, 2015) pp. 451– 455. [24] G. Zaslavsky and M. Edelman, Chaos: An Interdisci- plinary Journal of Nonlinear Science 10 , 135 (2000). [25] A. W olf, J. B. Swift, H. L. Swinney , and J. A. V astano, Ph ysica D: Nonlinear Phenomena 16 , 285 (1985). [26] G. L. Kotkin and V. G. Serbo, Col le ction of Pr oblems in Classic al Me chanics , V ol. 31 (Pergamon Press, 1971) pp. 179–182. [27] R. Mon tgomery , American Mathematical Mon thly 122 , 299 (2015). [28] T. Iw ai and H. Y amaok a, Journal of Ph ysics A: Mathe- matical and General 38 , 5709 (2005). [29] R. G. Littlejohn and M. Reinsch, Reviews of Mo dern Ph ysics 69 , 213 (1997). [30] M. Stomp, J. Huisman, L. J. Stal, and H. C. Matthijs, The ISME journal 1 , 271 (2007). [31] D. Bambusi, Univ ersita di Milano (2014). [32] D. Shep elyansky , Physical Review E 82 , 055202 (2010). [33] A. S. De Castro and M. F. Suga ya, Europ ean journal of ph ysics 14 , 259 (1993). [34] L.-S. Y ao, Nonlinear Analysis: Mo delling and Control 15 , 109 (2010).