{'fraction_non_alphanumeric': 0.03943355119825708, 'fraction_numerical': 0.014911643669813604, 'mean_word_length': 4.360887620036335, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bottlenecks to vibrational energy flow in OCS: Structures and mechanisms

Finding the causes for the nonstatistical vibrational energy relaxation in the planar carbonyl sulfide (OCS) molecule is a longstanding problem in chemical physics: Not only is the relaxation incomplete long past the predicted statistical relaxation time, but it also consists of a sequence of abrupt transitions between long-lived regions of localized energy modes. We report on the phase space bottlenecks responsible for this slow and uneven vibrational energy flow in this Hamiltonian system with three degrees of freedom. They belong to a particular class of two-dimensional invariant tori which are organized around elliptic periodic orbits. We relate the trapping and transition mechanisms with the linear stability of these structures.


I. INTRODUCTION
How does vibrational energy travel in molecules? Answering this question succinctly seems a hopeless task considering the complexity of interatomic interactions in a molecule. Yet even before scientists were burdened by this knowledge, the so-called statistical theories posited the answer: Vibrational energy travels "very fast" and distributes itself statistically among the vibrational modes of a molecule, assumed to resemble an assembly of coupled oscillators, well before a reaction takes place. Reaction rate theories based on these assumptions -known collectively as statistical or RRKM theories [1,2,3,4] -remain reliable working tools of the practicing chemist because they have been vindicated in an overwhelming number of chemical reactions.
However, numerical studies of Hamiltonian systems have provided solid evidence [5,6,7,8,9,10] that the approach to equilibrium usually proceeds more slowly than predicted by statistical theories [11,12] -and it is also nonuniform, showing intriguing fits and starts. In particular, for Hamiltonian systems with two degrees of freedom, the familiar picture of chaotic seas, rigid boundaries in terms of noble tori [13], leaky barriers in terms of cantori [14,15] has been well-established in the literature, and these structures are found to be the source of anomalous transport in such systems [16,17].
Beyond two degrees of freedom, the transport picture in terms of phase space structures is less clear. However, the phase space of higher-dimensional systems shows similar features such as the abundance of periodic orbits, and a mixture of chaotic and regular regions, the latter being characterized (under some hypothesis) by invariant tori of various dimensions. The KAM theorem [13] states that these structures are in general robust with respect to an increase of the perturbation or equivalently to an increase of energy. Understanding transport properties has to rely on these robust structures which are encountered by any typical trajectory. Roughly speaking, the presence of so many periodic orbits explains why generic trajectories, even when the system is strongly chaotic, display long intervals of near-regular behavior alternating with fits of chaos-a hallmark of anomalous diffusion.
The slow approach to equilibrium started to be acknowledged a little over fifty years ago with the investigation of the dynamics of coupled oscillators by Fermi, Pasta and Ulam who showed that the relaxation problem is far more complex than anticipated [5,6,7,8,9,10]. In chemical physics, anomalous diffusion was first implicated in the intramolecular vibrational energy relaxation of the carbonyl sulfide OCS molecule [11]. The numerical study of a classical Hamiltonian model of OCS shows very slow energy redistribution among the vibrational modes, even in the fully chaotic regime [11], disagreeing strongly with the fast timescales derived from traditional statistical theory. The understanding of the dynamics was successfully achieved for a collinear model of OCS which has two degrees of freedom [18,19,20,21,22]. However, severe technical difficulties [23,24,25] have prevented such a level of understanding beyond two degrees of freedom, and in particular, for the planar OCS model, in which the molecule is allowed to bend.
In this paper, we analyze the dynamics of a model for the planar OCS which is a Hamiltonian system with three strongly coupled degrees of freedom. The aim is to identify the relevant structures in phase space which are responsible for trappings and escapes, strongly influencing the transport properties (most prominently, the redistribution of intramolecular energy among the three modes). For example, rapid diffusion through phase space takes place through the so-called accelerator modes [26]. In contrast, sticky structures [27] like resonant islands or tori influence the dynamics by strongly slowing down the trajectories passing nearby. All these structures are responsible for anomalous diffusion and fractal kinetics in the system (for recent surveys, see Refs. [16,17] and references therein). Identifying these structures and the mechanisms behind trapping, escape and roaming is essential for understanding the transport properties of a given system. Given that there are many such structures in a realistic system, the only realistic hope for forming a generally valid picture of transport is to locate invariant structures which are responsible for the main changes in the transport properties.
The specific question we address is: What are the structures in the phase space of OCS acting as dynamical bottlenecks to the diffusion of chaotic trajectories? What are the structures allowing transitions to other parts of phase space? For three degree of freedom systems, these invariant structures can be invariant tori with dimensions zero (stagnation points), one (periodic orbits), two or three [28,29,30]. They can also include the stable and unstable manifolds of these objects [31]. How are invariant structures relevant in the phenomena of capture in chaotic systems? For planar OCS, we find that the bottlenecks and the transition mechanisms from trapped to hyperbolic behavior are provided by a particular class of two-dimensional tori and their unstable manifolds. These results were recently announced in a Letter [32].
The paper is organized as follows: In Sec. II, we briefly recall some basics of the Hamiltonian model for the planar rotationless OCS molecule. We also summarize the main results obtained on the dynamics of OCS relevant to the transport properties (both in the planar and collinear cases). In Sec. III A, we illustrate the transitions which occur in the neighborhood of periodic orbits using several representations: Time series, time-frequency analysis, and Poincaré sections. The striking common feature exhibited by many trajectories support the idea of some kind of universal transition mechanism. In Sec. III B, after summarizing our methodology, we investigate the neighboring phase space structures which strongly influence the dynamics of these trajectories.

A. The Hamiltonian
The dynamics of the planar model of carbonyl sulfide (OCS) can be described by a Hamiltonian model with three degrees of freedom with three strongly coupled, non-separable modes: There are two stretching modes and one bending mode. Each mode is represented by a coordinate and momentum pair, which we define as: R 1 = d(C, S) and P 1 for the CS stretching mode, R 2 = d(C, O) and P 2 for the CO stretching mode, and finally, α = ∠(OCS) and P α for the bending mode of the molecule (see Fig. 1). Hamiltonian model for the rotationless OCS molecule has been provided in [11,33] has the form where T is the kinetic energy and V is the potential energy. The kinetic energy is quadratic in the momenta and is provided as where µ i are the reduced masses. Based on available experimental data, the analytic model for the potential energy surface has been proposed in [33]. In summary, V is given by where V i (R i ) are Morse potentials for each of the three interatomic distances R 1 , R 2 and R 3 = d(S, O), and Here, R 0 i are the equilibrium interatomic distances, and R 3 is given by At equilibrium, the molecule is collinear, therefore R 0 3 = R 0 1 + R 0 2 . Also, the interaction potential V I assumes the Sorbie-Murrell form: where P (R 1 , R 2 , R 3 ) is a quartic polynomial in each of its variables: ij , c ijk , c (4) ijkl ) are provided in Ref. [11]. We display the equipotential surfaces of V (R 1 , R 2 , R 1 + R 2 ) of the collinear configuration in Fig. 1. The equations of motion can be derived from Hamiltonian (1) using the canonical Hamilton's equations.

B. Summary of prior results on the OCS dynamics
The classical models of both the collinear and the planar (rotationless) carbonyl sulfide OCS molecule have been studied in detail in Refs. [11,18,21,31,34,35,36].
The dynamics in the collinear configuration of OCS was first studied by Carter and Brumer [11]. They characterized the motion of this system at a number of energies, extending up to 20, 000 cm −1 (which amounts to E = 0.09 a.u.) A relaxation time, as defined in Refs. [37,38,39], was estimated at 0.17 pico-seconds. However, after integrating trajectories for 2.4 picoseconds, no relaxation to statistical equilibrium was observed. When this contradiction was investigated by integrating the equations for much longer times (up to 45 picoseconds), two distinct timescales for relaxation were found, the longer of which characterized energy redistribution that was incomplete even after 45 picoseconds [35]. Even on the picosecond time scale, sudden transitions between relatively long-lived regions of localized mode energies were observed. Since this collinear model has two degrees of freedom, Davis and Wagner [35] used Poincaré surfaces of section as a visualizing tool for phase space structures. These revealed that even at high energy (E = 0.09), the system has a "divided phase space", with coexisting regular and chaotic regions. They observed that trajectories can be trapped in restricted regions of phase space for many vibrational periods, after which they would suddenly move to other regions of phase space to repeat the pattern.
Progress came with the recognition that the thenrecent lobe dynamics [14,15] could help to explain non-statistical relaxation in two degree of freedom systems [18]. When the strength of the perturbation (or equivalently, the total energy) is increased, the twodimensional invariant tori of a Hamiltonian system with two degrees of freedom develop sets of "holes" with the systematics of Cantor sets. These holes, dubbed "cantori" [14], form leaky barriers which can act as bottlenecks to phase space transport. These bottlenecks are associated with broken tori with irrational frequency ratios, where those with "noble" number ratios being generically the very last to be destroyed by an increasing perturbation (the supporting argument being that these numbers are the most poorly approximated by rationals [40]). For OCS, their existence has been confirmed in Ref. [18] in a region between two resonances ω CO /ω CS = 3/1 and ω CO /ω CS = 5/2. The noblest irrational number between the rationals 5/2 and 3/1 is 2+γ, where γ = ( √ 5−1)/2 is the golden mean [14,40] and can be expressed as a continued fraction of an infinite sequence of ones, also written as . These results obtained from classical mechanical were confirmed using quantal wave packet calculation [41]. However, these successful results could not be extended to the planar OCS due to severe technical and computational difficulties [23,24]. Yet there were indications that this problem of intramolecular energy flow in higher dimensions is also related to the resonant and non-resonant structures [21,36]. In particular, the relevance of Arnold's web in the diffusion of trajectories was highlighted. Among their conclusions are that transport is most rapid along low order resonance zones; transport is slow (diffusive) along high order resonances; it was conjectured that pairwise noble frequency ratios play a role of inhibiting transport along resonance lines. The complexity of transport processes in the collinear OCS model, revealed in the early investigations, suggests that a look into phase space structures such as periodic orbits or invariant tori is needed for a better understanding of these processes. Even if the measure of such invariant structures embedded into a chaotic sea is typically zero, the "neighborhoods" of influence around them can have relatively large measure and their finite-time properties, as characterized by Lyapunov exponents, provide a quantitative picture of transport. The rationale goes as follows : An ensemble of trajectories, described by a density function, which is centered in a finite volume around a periodic orbit, will evolve in finite time following this periodic orbit, and spreading predominantly in the direction of unstable manifolds, exponentially in time with a rate equal to the local Lyapunov exponent. An orbit in the "neighborhood" of a periodic orbit, temporarily assumes or "shadows" the properties of this periodic orbit as a general consequence of dynamical continuity [42]. This temporary influence of periodic orbits can also be viewed as instantaneous time-periodic forcing, exerted by a periodic orbit. It is expected that, in general, typical trajectories are trapped for longer times in the neighborhoods of linearly stable orbits. In what follows, we draw a dynamical picture of transport in OCS based on the determination of invariant structures in phase space and their linear stability properties.

Density of periodic orbits
A generic feature of Hamiltonian dynamics is the abundance of periodic orbits in phase space. Figure 2 represents the averaged density in the configuration space (R 1 , R 2 ) of periodic orbit points on the Poincaré section Σ (defined in Sec. III A 3) for planar OCS. A closer inspection of this figure shows that the most prominent regions of stability surround short periodic orbits with elliptic linear stability. A typical trajectory passes through this maze of periodic orbits, being trapped for some time according to local stability properties. The aim of this manuscript is to understand how a typical trajectory can be trapped and released locally around a given periodic orbit. In what follows we analyze the transport properties in the neighborhood of an elliptic periodic orbit, like for instance O a , as shown by a circle in Fig. 3.

Time-frequency analysis and stroboscopic mapping
To examine the temporal features of trajectories we use time-frequency analysis [43].
In what follows, x designates a point in phase space, i.e. x = (R 1 , R 2 , α, P 1 , P 2 , P α ). A finite segment of a trajectory can be represented by a sequence of phase space coordinates, {x n } n=1,...,N , x n = x(t n ), visited by a trajectory at times t n . For a stroboscopic map, we take snapshots with a fixed time increment, t n+1 = t n + ∆. It is natural to scale the time increment ∆ by the period T p of the organizing periodic orbit. We select ∆ = T p /4. The time series of selected orbits are displayed in the bottom panels of Figs. 4 and 6.
We study the instantaneous frequencies using wavelet decomposition. As described in Refs. [43,44,45], the time-frequency analysis is based on a continuous wavelet transform of an observable f (t) We choose the mother wavelet ψ, in the Morlet-Grossman form: ψ(t) = e ιηt e −t 2 /2σ 2 /(σ 2 π) 1/4 , with adjustable parameters η and σ. The time-frequency representation is obtained via a relation between the scale s and the frequency ξ = η/s. We consider the normalized scalogram which can be interpreted as the energy density in the time-frequency plane. The ridges of P W can be interpreted as instantaneous frequencies, or more rigorously,   the set of frequencies for a given time interval. In this section, two typical trajectories (whose initial conditions are specified in Tab. I) which are initially close to elementary organizing periodic orbits, are represented in Figs. 4 and 6 where the signal f (t) is chosen to be P 1 (t) or P 2 (t). It should be noticed that other choices of ob-  servables f (t) lead to the same qualitative features as the ones presented here, and in addition, these features are common to a wide set of other trajectories in the same neighborhood.
Time-frequency analysis shows that each of these trajectories displays qualitatively distinct regions : some with approximately constant ridges in time, and others   with multiple and short ridges. These two regimes are clearly marked with transition intervals (highlighted by shaded bands in the figures). In Fig. 4 we observe a clear, sharp transition stage between trapped behavior (around O a ) and roaming behavior throughout a large portion of phase space. After some time spent around the periodic orbit, the trajectory seems to find an exit channel through a bottleneck. Generically, any trajectory experiences multiple events of capture and escape (like the one in Fig. 6). We have found that escape to the chaotic region proceeds in two stages, characterized by two different rates of escape. The transition interval is characterized by t trap and t esc . The first ("slow") stage, 0 ≤ t ≤ t trap , and the second ("fast") stage t trap ≤ t ≤ t esc . The precise transition points located at t trap and t esc may vary in different situations, but typically t esc − t trap ≪ t trap . We notice that these trapping and transition stages although visible, were not as clearly apparent on the time series as on the time-frequency plots.
In order to identify the phase space regions visited during the trapping and escape stages, we complement the time-frequency analysis by projections of segments of the time series in a two-dimensional plane given by two coordinates, e.g., (R 1 , R 2 ), (R 1 , P 1 ) or (R 2 , P 2 ). In Fig. 7, two segments of the trajectories of Fig. 4 are represented (left and middle panels), one corresponding to the trapped stage (to the left of the shaded band), and the other one to the trajectory after the escape process (to the right of the shaded band). It is shown that the trajectory is trapped into a small L-shaped region around the periodic orbit O a , and that after the transition point, the trajectory has access to a larger part of phase space with an apparent size of the order of the entire accessible region. The same observation follows for the trajectory of Fig. 6, the stroboscopic plot of which appears in Fig. 8. By drawing tubes around it, we notice that the trajectory in both trapped segments sticks to particular regions around different periodic orbits.

Poincaré sections
We use Poincaré sections as another way to visualize multidimensional trajectories. Given that this Hamiltonian system has three degrees of freedom, the Poincaré section is four dimensional. We show below how twodimensional projections of these sections can be used to gain insight into the dynamics (although this information is displayed less clearly for this system than for a system with two degrees of freedom). Given some scalar function U (x) of the phase space variables, we define this section Σ to be the set of points x of a trajectory such that U (x) = 0 , withẋ · ∂U/∂x > 0. From two consecutive points x n = x(t n ) and x n+1 = x(t n +∆(x n , t n )) on the Poincaré section, we define a Poincaré map F Σ , In what follows, we have used the surface Σ defined by The argument for choosing this surface goes as follows: We study an energy range where the time series of the bending mode (α, P α ) oscillate around an instantaneous mean value α . Between each oscillation there is a turning point where momentum P α vanishes. The only case which is not captured by the Poincaré section is when the bending mode is "frozen" to α = π and P α ≡ 0, which corresponds to the collinear OCS.
We choose a four-dimensional parametrization of the surface of section Σ which consists of R 1 , P 1 , R 2 and P 2 . Setting P α = 0 in Eq. (6), the equation with a constraintṖ α > 0 is to be solved for α(R 1 , P 1 , R 2 , P 2 ; E) numerically. There are two merits in using Poincaré sections: First, representing projections as a set of planar plots of canonically conjugated variables helps in perceiving the symplectic symmetry of structures. Second, the section manifold is one dimension smaller than the energy manifold, and so are the maps of all the invariant structures. For instance, periodic orbits correspond to a finite set of points {x n } n=1,...,N on Σ, and the dynamics visits these points in a cyclic manner, i.e. F Σ (x n ) = x n+1 for n = 1, . . . , N − 1, with F Σ (x N ) = x 1 . Similarly two-dimensional tori correspond to closed curves on Σ. In Figs. 9 and 10, Poincaré sections, projected on the planes (R 1 , P 1 ) and (R 2 , P 2 ), are drawn for the two trajectories considered in Figs. 4 and 6. These Poincaré sections clearly show distinct one-dimensional curves (clearly visible in the insets of Figs. 9 and 10) in the transition stages (shaded bands on Figs. 4 and 6). The tubes which we identify as two-dimensional invariant tori in phase space [32], represented on Figs. 7 and 8, correspond to these one-dimensional curves (or more generally a set of one-dimensional curves) on the Poincaré sections.
In the trapping stages (around specific periodic orbits) like the ones in Figs. 7 and 8 (left panels), the rotation numbers are obtained from the frequency map analysis [46] on the surface of section Σ. Dimensionless ratios of frequencies arise naturally in the Poincaré map F Σ , and ratios of frequencies are called rotation numbers. The trapping stage can be characterized by a single rotation number (and its harmonics), implying that a twodimensional torus is the relevant invariant structure in the trapping process. In the following, we determine such structures and highlight the family of two-dimensional tori which are relevant for the transport picture in this system.

A summary of the methodology
It is well established that invariant structures in phase space play an important role in the transport properties associated with Hamiltonian systems with two degrees of freedom [13]. In particular, the role of periodic orbits has been singled out in many experiments [47]. Even if some aspects of this dynamical picture can be extended to systems with a larger number of degrees of freedom, it remains to address the role of invariant structures which are not present in systems with one and two degrees of freedom, but are specific to three and more degrees of freedom. In three degree of freedom systems this new type of invariant structures takes the form of two-dimensional invariant tori. Observations described in Sec. III A indicate that such tori close to elliptic periodic orbits play an important role. To have a qualitative description of dynamics near a periodic orbit, we consider a fixed point x 0 on the surface of section Σ, i.e. F Σ (x 0 ) = x 0 , corresponding to a point of a periodic orbit. Near x 0 , the Poincaré map F Σ can be expanded into a linear part and a remainder: where DF Σ (x 0 ) is the matrix of first order derivatives of the Poincaré map, constrained to the surface of section and evaluated at x 0 . All higher order terms in x − x 0  Fig. 4 (for E = 0.09) on the (R1, P1)-plane (left panel) and on the (R2, P2)-plane (right panel). It is apparent (in the insets) that before escaping to the external region, the trajectory is stuck in the neighborhood of a well-localized structure. The five points (in bold) are the intersections with Σ of a partially hyperbolic resonant periodic orbit which is responsible for the escape to the chaotic region through its unstable manifold. Two "tentacles" starting at the upper and lower parts of the structure are marked with broken lines connecting crosses to clarify what happens during the escape stage (shaded band in Fig. 4). The resonant 2:5 periodic orbit has elliptic-hyperbolic stability with λ = 0.113231998 per return (or λ = 0.566159992 for the entire orbit) and a rotation number ω/π = 0.35684077865194591 per entire orbit. are collected in R(x−x 0 ). Finite-time dynamics near the fixed point x 0 are determined by the properties of the matrix DF Σ (x 0 ). Assuming that linearized approximation is effective, and discarding the remainder term from further discussions (the fully nonlinear problem with large R is solved using the methodology outlined in Appendix B), we consider a closed curve γ(s) on the Poincaré section Σ defined on a torus s ∈ T 1 , and consider the dynamics of x(s) = x 0 + ǫγ(s) given by If F Σ has at least one pair of eigenvalues in the form Λ = exp (±ιω), it is possible to find a γ(s) such that DF Σ γ(s) = γ(s + ω ǫ ) and |ω − ω ǫ | = o(ǫ). Therefore the equation has a family of solutions, parametrized by the rotation number ω. Equation (9) defines a torus as a loop on the surface of section Σ with rotation number ω. Even if DF Σ (x 0 ) has two pairs of eigenvalues of the form exp (±ιω i ) such an invariant loop close to x 0 can be found. More details on the determination of twodimensional invariant tori are given in Appendix B.

Invariant tori and their bifurcations: Bottlenecks
In the cases discussed in Sec. III A, trajectories undergo a transition (after a trapping stage) in the vicinity of a nonresonant elliptic periodic point x 0 , whether it is associated with O a or O b . For each of these periodic points, the matrix DF Σ (x 0 ) has eigenvalues exp (±ιω 1 ), exp (±ιω 2 ) (numerical values are given in Tab. I.) Processes associated with the escape from the trapping stage can be better understood by analyzing the tangent space of the elliptic periodic orbit O a that locally has the structure of a direct product (center + center) T 1 ×I 1 ×T 1 ×I 2 , with the periodic orbit at the origin. The elements of the two intervals I i ⊂ R are rotation numbers ω i , which are not unique in general: The choice is fixed by requiring lim ǫ→0 ω i = ω 0 i , where ǫ is a measure of the "diameter" of the torus and ω 0 i are stability angles of the organizing periodic orbit. The Poincaré map F Σ induces rotations on T 1 , r ω1 × 1 × r ω2 × 1, where r ω is a rotation on T 1 with the rotation number ω. Partial (or complete) resonances are determined by one (or two) resonance conditions nω 1 + mω 2 + k = 0, where (n, m, k) are integers such that |n| + |m| + |k| > 0. The most striking trapping effects are observed for partial resonances of the type They are two-dimensional manifolds (locally), and can be foliated by one-dimensional invariant closed curves, called hereafter "loops." We propose to investigate a resonance manifold by mapping out dynamical invariants that form its backbone structure. Choosing either of the two situations, a resonance channel has been constructed by finding the two-dimensional invariant tori for ω i ∈ I i . At a small distance from the periodic orbit we use information obtained from the linear normal form DF Σ (x 0 ). Once an initial loop is found, we follow the progress as the rotation number ω is varied continuously monitoring their stability properties. Local normal stability of each family of tori can be represented by plotting the maximal Lyapunov exponent λ by solving the generalized eigenvalue problem [See Appendix B and Eq. (B3)], versus the rotation number ω. Such a plot for a family of two-dimensional tori, originating from O a , is shown in Fig. 11. From Fig. 11, we obtain a transition point at ω/π ≈ 0.481 in the form of a bifurcation of an invariant torus. Projections of two-dimensional invariant tori in the transition regions are shown in Figs. 7 and 8, while corresponding loops in the surface of section Σ are shown in Figs. 9 and 10.
The transition stage (see Figs. 9 and 10) indicates an exponential divergence resulting from an escape along the unstable branch of a hyperbolic manifold. The maximal Lyapunov exponent of the segment of a trajectory in the capture stage can be estimated by observing the duration of capture, and the per-return Lyapunov exponent can be estimated as λ ≃ 1/N , where N is the number of returns to the surface of section Σ before the escape. In Fig. 9 the proximity of the torus to a 2:5 resonance zone suggests the influence of a periodic orbit with or 2:5 winding number ratio (in hollow circles). For trajectory close to orbit O a we have N ≈ 150, yielding a typical value of λ ≃ 0.007. This value is inconsistent with the Lyapunov exponent of the nearby resonant periodic orbit (which has a Lyapunov exponent of λ = 0.113 per return to the surface of section Σ), indicating that other structures than periodic orbits are important in describing the capture processes. Unstable two-dimensional tori are indeed better candidates for the escape scenario : An estimate of the Lyapunov exponent in the escape stage is consistent with the scenario of escape along unstable manifolds of the resonant orbit. The local rate of transition at the onset is estimated by the largest Lyapunov exponent in the family. In the case shown in Fig. 11 it is close to λ = 0.06. The full picture of dynamics is complicated by existence of a family of tori with varying (and smaller) Lyapunov exponents.
In Fig. 12 we represent the two families of twodimensional tori considered in Fig. 11 (blue and black  curves). First, the organizing periodic orbits O a (center of Fig. 12) and the resonance 2:5 (exterior spheres) are located. The projections of the two families of loops in Σ are plotted in the three dimensional space (R 1 , P 1 , R 2 ). Meridians of the surfaces are invariant under the Poincaré map (i.e. they are invariant loops γ). The first family of tori (blue curve in Fig. 11) starts from the central periodic orbit O a and continues outwards as the rotation number decreases from the value of ω (Oa) /π = 0.49 (see Table I and the blue curve in Fig. 11). The first loops of this family have zero Lyapunov exponent (the ones with ω/π between 0.49 and 0.481). At the bifurcation point (ω/π = 0.481), the second family (black curve on Fig. 11) branches off of the first one and continues nor-mally (with zero Lyapunov exponent). The continuation of the first branch of tori is now normally hyperbolic (see Fig. 11) from ω/π between 0.481 to 0.432, while the new branch of frequency halved loops is at first elliptically stable. The bifurcation at ω/π = 0.481 is a frequencyhalving, since the emerging loop winds around the original one twice, or in other words, has half the rotation number. This process is very general and we expect it to occur in the vicinity of any periodic orbit with several elliptic stability degrees of freedom. The family of tori has singularities at some specific rotation numbers, but the manifold can typically be continued across them, and therefore seems to be robust. The behaviour of the second branch of this family of tori as it approaches the rational rotation number ω/π = 2/5 was investigated. A nontrivial foliation of invariant loops in the vicinity of a 2:5 periodic orbit is shown in Fig. 13. Allowed by dimensional analysis, a possible scenario is that this family of tori is heteroclinic to the invariant manifolds of other invariant tori, related to periodic points in 2:5 resonance with O a . However, a picture of interconnected families of tori, permeating bulk of the entire phase space is yet to emerge.
From the numerical simulations of a large assembly of trajectories, the following assumptions emerge : 1) the lowest order k = 1 resonance controls the rates of transition from regular to chaotic dynamics, 2) the k = 1 resonance is a manifold that has a two-dimensional "backbone" manifold, in analogy with resonance manifolds of integrable Hamiltonian systems, and 3) regular-tochaotic transition occurs at the point where there is a transition in the normal stability of this manifold. From these assumptions, the typical scenario for escape after trapping by a weakly hyperbolic family of tori, is the following one: First the trajectory evolves in a regular region until it finds an exit channel (the transition stage) in the form of a manifold of normally hyperbolic invariant two-dimensional tori, and follows along a manifold becoming more chaotic progressively, as it visits invariants with larger hyperbolicity (Lyapunov exponent). Eventually it is escapes to a strongly chaotic region using the unstable manifold of a hyperbolic periodic orbit with a large Lyapunov exponent.

IV. CONCLUSIONS
In contrast to collinear OCS where the phase space is roughly divided into islands and chaotic seas, the phase space of planar OCS exhibits a complex ocean with currents, reefs and shoals which slow down the progress toward energy equilibration. In this article, we have identified these structures and their linear stability properties. Principal among them are two-dimensional invariant tori which occur in families and can be parametrized by their rotation numbers. These structures are organized around periodic orbits which provide the backbone to the dynamics. By trapping trajectories temporarily, they act The first branch (nearly horizontal and corresponding to the blue curve in Fig. 11) emerges close to the fixed point Oa (in the center of the figure). The second branch (nearly vertical and corresponding to the black curve in Fig. 11) emerges at the bifurcation point of the first branch (with ω/π = 0.481422634). We also represent the points (in black) of the 2:5 resonant periodic orbit which obstructs the continuation of the second branch. Note that this figure links as bottlenecks to the exploration of larger parts of phase space. Our work also makes explicit the mechanisms by which trajectories are trapped and by which they escape from the trap. FIG. 13: Approach to a rational rotation number 2:5 in the second branch of a family of tori of periodic orbit Oa, E = 0.09 (see also Fig. 11.) Rotation numbers, from right to left are 2ω/π = 0.4148,0.4110,0.4104. Nontrivial foliation around manifolds of a resonant 2:5 periodic orbit can be seen. Note that this figure is related to frequency-halved loops in the second branch of Fig. 11, therefore we have doubled the original rotation numbers. classification) symmetry group C 2v which acts on intrinsic coordinates P 1 , P 2 , α and P α , while R 1 and R 2 are left invariant. Elements of C 2v are identity e, reflection σ 1 , reflection σ 2 , and inversion i, defined as : e(P 1 , P 2 , α, P α ) = (P 1 , P 2 , α, P α ), σ 1 (P 1 , P 2 , α, P α ) = (P 1 , P 2 , 2π − α, −P α ), σ 2 (P 1 , P 2 , α, P α ) = (−P 1 , −P 2 , 2π − α, P α ), i(P 1 , P 2 , α, P α ) = (−P 1 , −P 2 , α, −P α ).
This discrete symmetry is useful for the method of surface of section, because it allows to relate points x in phase space with P α = 0 andṖ α < 0, not on the surface Σ, with points σ 1 (x) which are on the surface Σ.
In addition to exact discrete symmetries discussed above, the specific form of potential energy (2) induces an approximate R 1 -R 2 reflection symmetry as seen in Fig. 1. Equation (2) can be written in the form of This partition quantifies the approximate symmetry. The non-vanishing parameters δD, D 3 and A measure the deviation from the exactly symmetry. Considering the linear transformations of the coordinates R 1 and R 2 given by L(R 1 , R 2 ) = (a 1 R 2 + b 1 , a 2 R 1 + b 2 ), the symmetry line is obtained by requiring that The solution is obtained in terms of parameters β 1 and β 2 , and in particular β 1 /β 2 ≈ 0.9, and the symmetry is then given by the equation : In case of an exact symmetry, the symmetry line would be a natural boundary of the elementary cell of the dynamics. All orbits could be classified with respect to this symmetry as having a symmetric partner, or being self symmetric, as usually. When the symmetry is only approximate the cell boundary argument is no longer valid, but the orbits can still be classified in this way, in particular, with regards to their degeneracy.

APPENDIX B: METHODOLOGY: DETERMINATION OF INVARIANT TORI AND THEIR LINEAR STABILITY PROPERTIES
We briefly summarize the method we used to compute two dimensional invariant tori of a Hamiltonian system. We have seen that this is equivalent to determining closed invariant curves (loops) of the Poincaré map on the chosen surface of section Σ. Furthermore we compute the linear stability properties of such objects. This method follows the one described in Ref. [48].

Determination of invariant tori
In order to determine two-dimensional tori, we use the fact that the type of internal dynamics on T 1 is likely to be a rotation. We assume that the Poincaré map F Σ has an invariant curve with an irrational rotation number ω, and that there exists a map (at least continuous) x : T 1 → Σ such that Denjoy's theorem [49] states that such a rotation number ω can be defined. Let C(T 1 , Σ) be the space of continuous functions from T 1 in Σ, and let us define the linear map T ω : C(T 1 , Σ) → C(T 1 , Σ) as the translation by ω, i.e. (T ω x)(θ) = x(θ + ω). We define F : It is clear that zeros of F in C(T 1 , Σ) correspond to (continuous) invariant curves of rotation number ω. The determination of two-dimensional invariant tori follows from the search of zeros of this functional. First we expand x(θ) in a Fourier series with real coefficients, x(θ) = a 0 2 + k>0 (a k cos πkθ + b k sin πkθ) , where a k , b k ∈ R n for k ∈ N (n being the dimension of the flow) and x(θ) is a periodic function with period 2, i.e. x(θ + 2) = x(θ). We truncate these series at a fixed value of N , and determine an approximation to the 2N +1 unknown coefficients a 0 , a k , and b k for 1 ≤ k ≤ N . We construct the discretized version of Eq. (9) by considering a mesh of 2N + 1 points on T 1 : where we notice that for numerical stability reasons, the length of T 1 is taken as 2. Given the Fourier coefficients a k , b k , the coordinates x(θ j ) are expressed as linear functions of the coefficients a k , b k , i.e.
x(θ j ) ≡ φ({a k }, {b k }, j), given by Eq. (B2). Accordingly, F Σ (x(θ j )) and Eq. (9) can be considered as functions of the coefficients a k , b k : for 0 ≤ j ≤ 2N and where i(ω) = (2N + 1)ω/2. The coefficients a k , b k are the unknowns in the above equation. We solve F = 0 using a Newton's iterative algorithm. At each iteration, it provides the corrections δa k and δb k to be added to the a k and b k obtained from the previous iteration. We approximate (δa, δb) as a solution of the following equation: where a = (a 0 , a 1 , . . . , a N ) and b = (b 1 , . . . , b N ). The iteration a ′ = a + δa, b ′ = b + δb and ω ′ = ω + δω converges if the initial guess is close enough to the true solution. The above equation requires the inversion of the Jacobian of F j . From the previous definitions it is clear that if x(θ) is a Fourier series corresponding to an invariant curve then, for any ϕ ∈ T 1 , y(θ) ≡ x(θ + ϕ) is a different Fourier series corresponding to the same invariant curve as x(θ). This implies that the Jacobian of F j around the invariant curve has, at least, a onedimensional kernel. To solve this problem we use the Singular Value Decomposition. Even if Newton's algorithm has converged, we cannot claim with certainty that a smooth two-dimensional torus has been found. We have noticed that crude discretization can wash out the details of non-smooth curves. Sometimes doubling the number of points in the discretization turns a convergent case into a divergent one. In most cases the reliability of a solution is almost certain when testing the spectrum of the solution (and the norm of its eigenvectors weighted by the frequency, penalizing high harmonics): a smooth solution should contain a unit eigenvalue. This is why it is also important to monitor the linear stability properties of the curves we obtain numerically.

Linear stability properties
In addition to the determination of the location of the invariant tori, we compute their linear stability properties to obtain information on the dynamics in its (infinitesimal) neighborhood, i.e. eigenvalues and eigenvectors which give at first order an approximation to the invariant manifolds (stable, unstable and central) near the invariant curve. We consider the generalized eigenvalue problem which amounts to finding (Λ,ψ) such that DF Σ (x(θ))ψ(θ) = Λψ(θ + ω). (B3) The eigenvalues Λ have the following properties [48]: 1) Λ = 1 is an eigenvalue of Eq. (B3); the corresponding eigenvector is the derivative of the loop x, 2) if Λ is an eigenvalue of Eq. (B3), then Λ exp (2ιkπω) is also an eigenvalue for any k ∈ Z, 3) the closure of the set of eigenvalues of Eq. (B3) is a union of circles centered at the origin.
There are two unit eigenvalues in the spectrum of DF Σ . The symplectic symmetry implies that the tori are degenerate in the linear approximation. It implies the existence of a family of (smooth) two-dimensional tori. As it is usual, we expect that this family is discontinuous and the discontinuities are around rational rotation numbers ω = m/n. Numerically, once an invariant torus with a specific ω is found, we simply increment the frequency parameter ω → ω + δω and restart the search. In this way, we determine these families of two-dimensional tori parametrized by their frequency ω on the Poincaré section. More details on the algorithm are given in Ref. [50].