Unraveling the escape dynamics and the nature of the normally hyperbolic invariant manifolds in tidally limited star clusters
MMNRAS , 525–546 (2017) Preprint 15 November 2018 Compiled using MNRAS L A TEX style file v3.0
Unraveling the escape dynamics and the nature of thenormally hyperbolic invariant manifolds in tidally limitedstar clusters
Euaggelos E. Zotos (cid:63) and Christof Jung † Department of Physics, School of Science, Aristotle University of Thessaloniki, 541 24, Thessaloniki, Greece Instituto de Ciencias F´ısicas, Universidad Nacional Aut´onoma de M´exico Av. Universidad s/n, 62251 Cuernavaca, Mexico
Accepted 2016 September 28. Received 2016 September 28; in original form 2016 August 30
ABSTRACT
The escape mechanism of orbits in a star cluster rotating around its parent galaxy in acircular orbit is investigated. A three degrees of freedom model is used for describingthe dynamical properties of the Hamiltonian system. The gravitational field of thestar cluster is represented by a smooth and spherically symmetric Plummer potential.We distinguish between ordered and chaotic orbits as well as between trapped andescaping orbits, considering only unbounded motion for several energy levels. TheSmaller Alignment Index (SALI) method is used for determining the regular or chaoticnature of the orbits. The basins of escape are located and they are also correlatedwith the corresponding escape time of the orbits. Areas of bounded regular or chaoticmotion and basins of escape were found to coexist in the ( x, z ) plane. The propertiesof the normally hyperbolic invariant manifolds (NHIMs), located in the vicinity ofthe index-1 Lagrange points L and L , are also explored. These manifolds are ofparamount importance as they control the flow of stars over the saddle points, whilethey also trigger the formation of tidal tails observed in star clusters. Bifurcationdiagrams of the Lyapunov periodic orbits as well as restrictions of the Poincar´e mapto the NHIMs are deployed for elucidating the dynamics in the neighbourhood of thesaddle points. The extended tidal tails, or tidal arms, formed by stars with low velocitywhich escape through the Lagrange points are monitored. The numerical results of thiswork are also compared with previous related work. Key words: stellar dynamics – galaxies: star clusters
Stars are born and developed from large clouds of moleculargas. These molecular clouds usually contain hundreds, oreven thousands, of solar masses of material, so the starsform in groups of clusters. After the remnant gas is heatedand blown away, the stars group together by their mutualgravitational attraction thus forming open star clusters.Dynamical interactions between the stars are the mainreason why star clusters have the tendency to gradually loosemass until they completely dissolved. A substantial amountof stars of the clusters become members of the disc andthe halo populations of the parent galaxy, as the clustersdissolve over time. The dissolution mechanism of star clus-ters is a very active as well as challenging field of research.Theoretical aspects of this intricate subject, like the evapo- (cid:63)
E-mail: [email protected] † E-mail:jung@fis.unam.mx ration process because of the movement of stars above theescape velocity and the time-scale of relaxation which deter-mines the rate of the dynamical evolution of a star clusterwere the first which have been investigated analytically (e.g.,Ambartsumian 1938; Spitzer 1940). Later on, the dynami-cal properties of stars which have been scattered above thecritical escape energy were examined (King 1959), while ayear later the importance of the close encounters betweenstars for the rate of mass loss of a star cluster was studied(H´enon 1960).According to Lada & Lada (2003) a considerableamount of the stars of a galaxy are actually born in starclusters. In particular, some recent studies shed some lighton the mechanism of how stars are in fact born and developstar clusters (e.g., Bressert et al. 2010; Kruijssen 2012). In-evitably, all star clusters dissolve over time due to two mainreasons: (i) the two-body process (encounters) that affectsall the members of a cluster thus forcing them to obtainescape velocities and (ii) the strong tidal forces due to the c (cid:13) a r X i v : . [ a s t r o - ph . GA ] F e b E. E. Zotos & Ch. Jung external gravitational attraction of the parent galaxy. Moreprecisely, the tidal field of the parent galaxy can significantlyaffect and therefore disturb the orbital behaviour of the starcluster itself (Ross et al. 1997).Usually we use the term “tidally limited” star clusterwhich means, roughly speaking, that the tidal field imposes aboundary, outside of which the mass density vanishes, whilethe star cluster is being captured within the tidal limit. Dur-ing the energy exchange between the stars, some of themreach the escape velocity and become runaway stars. Theescape of stars form a tidally limited cluster is a two-stageprocess. In the first stage, the stars are scattered into theescaping phase space by two-body encounters, while at thesecond stage they escape through the exit channels in theopen equipotential surface. The time required for a star tocomplete the first stage strongly depends on the relaxationtime, while on the other hand, the time needed for a star tocomplete stage two is mainly related to the energy the starhas. In Fukushige & Heggie (2000) an expression regardingthe escape time of stars in the second stage which have justfinished stage one had been derived, while in Baumgardt(2001) the corresponding results had been exploited in or-der to address the important issue of the dissolution time ofstar clusters. Furthermore, some additional interesting as-pects, like the velocity of escaping stars (e.g., Kollmeier &Gould 2007; Silva & Napiwotzki 2011; Li et al. 2012), ortheir population type (Royer 1997) have also been studied.Complicated stellar formations, such as tidal tails ortidal arms, are formed by stars which escape from the clus-ter and they are captured by the strong gravitational fieldof the parent galaxy. These interesting stellar formationshave been observed in the Milky Way (e.g., Grillmair etal. 1995; Kharchenko et al. 1997; Lehmann & Scholz 1997;Odenkirchen et al. 2001; Rockosi et al. 2002; Lee et al. 2003;Belokurov et al. 2006) and they have also been modelledand explored (e.g., Johnston et al. 1999; Yim & Lee 2002;Dehnen et al. 2004; Koch et al. 2004; Capuzzo Dolcetta etal. 2005; di Mateo et al. 2005; Lee et al. 2006; Choi et al.2007; Fellhauer et al. 2007; Montuori et al. 2007; K¨upperet al. 2008; Just et al. 2009; K¨upper et al. 2010). At thispoint we should emphasize, that all the above-mentionedreferences on tidal tails and tidal arms are exemplary ratherthan exhaustive.In Ernst et al. (2008) (hereafter Paper I) the escape dy-namics of a 2-dof model of a star cluster embedded in thetidal field of a parent galaxy was investigated, while in Zotos(2015a) (hereafter Paper II) the exploration was expandedinto three dimensions by using a 3-dof model. In both papersthe orbital dynamics and the corresponding basins of escapewere determined and obtained by conducting a thorough andsystematic orbit classification in the available phase space.In Zotos (2015b) the escape dynamics of a tidally limitedstar cluster was compared by using different types of spher-ically symmetric gravitational potentials for describing theproperties of the star cluster.If the energy is above, but close, to the saddle energyof the relevant saddle of the effective potential then the es-cape process is directed and channeled by the stable andunstable manifolds of some normally hyperbolic invariantmanifolds (NHIMs) sitting over this saddle. On this basis,a central part of the work of this article will consist in thestudy of the NHIMs. The Lyapunov orbits (Lyapunov 1907) are important parts of the NHIMs. Therefore, also the devel-opment scenario of the Lyapunov orbits will be investigatedin detail. We apply recently developed methods (e.g., Gon-zalez et al. 2014; Gonzalez & Jung 2015) for the numericalstudy of NHIMs and of the Poincar´e map restricted to theNHIM. This restricted map acts on a 2-dimensional domain,can be represented by 2-dimensional graphics and thereforeit is an ideal tool to present graphically the developmentscenario of codimension 2 NHIMs in 3-dof systems. In a re-cent paper (Jung & Zotos 2016) (hereafter Paper III) thedevelopment scenario of codimension 2 NHIMs in a 3-dofmodel of a barred galaxy has been numerically investigated.The structure of the present paper is as follows: In Sec-tion 2 we briefly describe the main properties of the tidallylimited star cluster model. In the following Section, we nu-merically investigate the escape dynamics of stars. Section4 contains a detailed description of the dynamics in theneighbourhood of the saddle points and in particular of theNHIMs sitting near these points. In Section 5 we link the in-variant manifolds with the tidal tail structures observed instar clusters. Our paper continues with Section 6, where themain conclusions of our work are presented. The paper endswith two appendices, in the first one we explain the connec-tion between the total potential and the effective potential,while in the second one we illustrate the scaling behaviournear the stable manifolds of NHIMs.
In order to describe the dynamical properties of a starcluster we need to use the so-called “tidal approximation”model. According to this theory we assume that the starcluster revolves around the center of a parent galaxy on acircular orbit and with constant angular velocity Ω. Usingthis assumption we are allowed to apply the epicyclic ap-proximation for determining the tidal forces which act onall stars of the cluster. For this purpose, the most appropri-ate system of coordinates is an accelerated rotating frame ofreference (Chandrasekhar 1942). In this system the center ofthe parent galaxy as well as the star cluster itself are at rest.Moreover, the origin of the coordinates is located at the cen-ter of the star cluster. This assures that the position of thegalactic center is P ( x, y, z ) = ( − R g , , R g is the radius of the circular orbit of the star cluster.Taking into account that the star cluster is sphericallysymmetric the best choice for describing its properties isby using the simple, yet self-consistent Plummer potential(more details about the isolated Plummer potential can befound in the Appendix A of Paper I)Φ cl ( x, y, z ) = − GM cl (cid:112) x + y + z + r , (1)where G is the gravitational constant, while M cl is the totalmass of the star cluster. The Plummer radius is r Pl = 0 . W = 4, where W isa parameter related to the mass density of the cluster (seee.g., King 1966).Considering the relationships connecting the tidal forceswith the epicyclic frequency κ as well as with the vertical MNRAS465
In order to describe the dynamical properties of a starcluster we need to use the so-called “tidal approximation”model. According to this theory we assume that the starcluster revolves around the center of a parent galaxy on acircular orbit and with constant angular velocity Ω. Usingthis assumption we are allowed to apply the epicyclic ap-proximation for determining the tidal forces which act onall stars of the cluster. For this purpose, the most appropri-ate system of coordinates is an accelerated rotating frame ofreference (Chandrasekhar 1942). In this system the center ofthe parent galaxy as well as the star cluster itself are at rest.Moreover, the origin of the coordinates is located at the cen-ter of the star cluster. This assures that the position of thegalactic center is P ( x, y, z ) = ( − R g , , R g is the radius of the circular orbit of the star cluster.Taking into account that the star cluster is sphericallysymmetric the best choice for describing its properties isby using the simple, yet self-consistent Plummer potential(more details about the isolated Plummer potential can befound in the Appendix A of Paper I)Φ cl ( x, y, z ) = − GM cl (cid:112) x + y + z + r , (1)where G is the gravitational constant, while M cl is the totalmass of the star cluster. The Plummer radius is r Pl = 0 . W = 4, where W isa parameter related to the mass density of the cluster (seee.g., King 1966).Considering the relationships connecting the tidal forceswith the epicyclic frequency κ as well as with the vertical MNRAS465 , 525–546 (2017) scape dynamics and NHIMs in star clusters frequency ν (see for more details Binney & Tremaine 2008)and also with the rotation of the star cluster we concludethat the total potential isΦ t ( x, y, z ) = Φ cl ( x, y, z ) + 12 (cid:0) κ − (cid:1) x + 12 ν z + Ω (cid:0) x + y (cid:1) . (2)Then the corresponding equations of motion read˙ x = p x + Ω y, ˙ y = p y − Ω x, ˙ z = p z , ˙ p x = − ∂ Φ t ∂x + Ω p y , ˙ p y = − ∂ Φ t ∂y − Ω p x , ˙ p z = − ∂ Φ t ∂z , (3)where, as usual, the dot indicates the derivative with respectto the time, while p x , p y and p z are the canonical momentaper unit mass, conjugate to x , y and z , respectively.For the computation of standard chaos indicators (theSALI in this case) we need to follow the time-evolution ofdeviation vectors w i , i = 1 ,
2. Therefore the required set ofvariational equations is˙( δx ) = δp x + Ω δy, ˙( δy ) = δp y − Ω δx, ˙( δz ) = δp z , ( ˙ δp x ) = − ∂ Φ t ∂x δx − ∂ Φ t ∂x∂y δy − ∂ Φ t ∂x∂z δz + Ω δp y , ( ˙ δp y ) = − ∂ Φ t ∂y∂x δx − ∂ Φ t ∂y δy − ∂ Φ t ∂y∂z δz − Ω δp x , ( ˙ δp z ) = − ∂ Φ t ∂z∂x δx − ∂ Φ t ∂z∂y δy − ∂ Φ t ∂z δz. (4)The equations of motion (3) admit the following iso-lating energy integral, which governs the motion of a testparticle (star) with a unit mass ( m = 1) H = 12 (cid:0) p x + p y + p z (cid:1) + Φ t ( x, y, z ) − Ω L z = E, (5)where E is the numerical value of the energy integral whichis conserved, while L z = xp y − yp x is the angular momen-tum along the z direction. At this point, it should be em-phasized that if the star cluster follows another type of orbit(i.e., elliptical instead of circular) around the parent galaxy,all numerical calculations become much more difficult to beperformed because no integral of motion, comparable to theenergy integral, is known.The whole system and the equations of motion of Eq.(3) have two discrete symmetries. First, they are invariantunder the transformation z → − z , p z → − p z , i.e. under a We would like to clarify that in Paper I the effective potentialΦ eff was presented, while now the total potential Φ t = Φ eff + (cid:0) x + y (cid:1) is used. At the end of the paper, in the Appendix A,we shall explain in detail how these two different concepts andtheir different functional forms are related. reflection in the symmetry plane z = 0. Second, the more in-teresting symmetry is a simultaneous inversion of all 4 hori-zontal coordinates, i.e. the transformation x → − x , y → − y , p x → − p x , and p y → − p y . This transformation is equivalentto a rotation of the whole system around the z -axis by anangle π .This second symmetry is a consequence of the tidalapproximation. As Eq. (2) shows some quadratic approx-imation has been made for the gravitational potential ofthe parent galaxy and the approximated potential containseven terms only. This property is the source of the symme-try. Any inclusion of odd terms into the total potential ofEq. (2) would destroy this symmetry. It is also obvious thatany two symmetry related orbits have the same value of theenergy because the Hamiltonian of the system is invariantunder the symmetry operation. Some consequences of thissymmetry will be discussed in later sections and it will bereferred to as the inversion symmetry of the system.For all the involved dynamical quantities we use thedimensionless system of units introduced in Paper I accord-ing to which G = 1, Ω = 1 and M cl = 4Ω − κ = 2 . r t = 1.The characteristic frequencies Ω, κ and ν on the other hand,are related to the galactic gravitational potential and in theneighborhood of our Sun can be expressed as functions ofthe Oort’s constants (more details can be found in Binney& Tremaine 2008). Using the numerical values of the Oort’sconstants provided in Feast & Whitelock (1997), we obtainthe following values: κ / Ω (cid:39) . ν / Ω (cid:39) . r t = (cid:18) GM cl − κ (cid:19) / . (6)The effective potential Φ eff has two Lagrange sad-dle points L and L which are located at L ( x, y, z ) =( − r L , ,
0) and L ( x, y, z ) = ( r L , , r L = 0 . . The numericalvalue of the effective potential at the Lagrange points yieldsa critical energy level E L = − . E = E L the equipotentialsurface encloses the critical volume, while for E > E L theequipotential surface is open and consequently stars are freeto escape from the cluster. In Fig. 1b we present a plot ofthe three-dimensional equipotential surface Φ eff ( x, y, z ) = E when E = − . > E L . We observe the two openings (exitchannels or throats) through which the test particles canleak out. In fact, we may say that these two channels aroundthe Lagrange points act as hoses thus connecting the inte-rior region ( − r L (cid:54) x (cid:54) r L ) of the cluster with the exteriorregion, or in other words the “outside world”. Exit channel1 (negative x -direction) indicates escape towards the galac-tic center, while channel 2 (positive x -direction) indicatesescape towards infinity.A double precision Bulirsch-Stoer FORTRAN 77 algo-rithm (e.g., Press 1992) was used for integrating backwards It is interesting to note that the Lagrange radius is almost equalto the tidal radius.MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung - - - - - - L L (a) (b) Figure 1. (a-left): The isoline contours of the effective potential Φ eff of the star cluster on the ( x, y )-plane for z = 0. The equipotentialcurve corresponding to the critical energy level E L is shown in red, while the position of the two Lagrange points L and L are indicatedby blue dots. (b-right): The three-dimensional equipotential surface of the effective potential Φ eff ( x, y, z ) = E when E = − .
23. The starsof the cluster can leak out through the exit channels 1 and 2 thus passing either through L or L , respectively. (For the interpretationof references to colour in this figure caption and the corresponding text, the reader is referred to the electronic version of the article.) and forwards the equations of motion (3) as well as the vari-ational equations (4). The time step of the numerical inte-gration was of the order of 10 − which is sufficient enoughfor the desired accuracy of our computations. Throughoutour calculations the numerical error regarding the conserva-tion of the energy integral of motion of Eq. (5) was smallerthan 10 − , although there were cases that the correspond-ing error was smaller than 10 − . All graphical illustrationspresented in this paper have been created using version 10.3of Mathematica (cid:114) (e.g., Wolfram 2003). In Paper II we investigated for the first time the escapedynamics of the three degrees of freedom system of the starcluster. In particular, we explored the orbital structure of theconfiguration ( x, y ) space as well as the phase ( x, ˙ x ) spacefor several initial values of the z coordinate. However it isnot possible to fully understand the escape process of starsby considering only some specific values of z .In this work in order to obtain a more complete viewregarding the escape mechanism of stars in tidally limitedstar clusters we shall conduct some additional numerical cal-culations. It would be very illuminating to unveil how theinitial value of z influences the escape process of the orbits.To achieve this we decided to define, for several energy lev-els E , dense uniform grids of 1024 × x, z ) plane inside the corresponding energetically al-lowed area. Following a typical approach, all orbits of the stars are launched with initial conditions inside the tidal ra-dius ( x + y + z (cid:54) r ). All three-dimensional orbits haveinitial conditions ( x , z ) with y = p x = p z = 0, whilethe initial value of p y is always obtained from the energy in-tegral (5) as p y = p y ( x , y , z , p x , p z , E ) (we choose thebranch of p y with positive orientation of intersection of theplane y = 0).Our task will be to distinguish between bounded andunbounded (escaping) motion. In Paper II we found a sub-stantial amount of trapped chaotic orbits, which do es-cape only after extremely long integration time. Therefore itwould be very useful to classify initial conditions of boundedorbits into two types: (i) non-escaping regular orbits and (ii)trapped chaotic orbits. Over the years several dynamical in-dicators have been developed for distinguishing between or-der and chaos. As in Paper II we choose to use the SmallerALingment Index (SALI) method (Skokos 2001) which hasbeen proved a very fast yet reliable tool. The nature of anorbit can be determined by the numerical value of SALI atthe end of the numerical integration. Being more precise,if SALI > − the orbit is ordered, while if SALI < − the orbit is chaotic. On the other hand, when the value ofSALI lies in the interval [10 − , − ] we have the case of asticky orbit and further numerical integration is needed tofully reveal the true character of the orbit.For the numerical integration of the initial conditionsof the orbits we set a maximum time of 10 time units.Our previous experience suggests that in most cases (energylevels) the majority of the orbits require considerable lesstime for finding one of the exit channels in the equipotential MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters surface and therefore escape (obviously, the numerical inte-gration is effectively stopped when an orbit escapes throughthe Lagrange points). Nevertheless, just for being sure thatall orbits have enough time to escape, we decided to usesuch a long integration time. Thus, any orbit which remainstrapped inside the tidal radius after an integration time of10 time units will be considered as a trapped (regular orchaotic) one.In the following we shall classify initial conditions oforbits into four categories: (i) trapped chaotic orbits, (ii)non-escaping regular, (iii) escaping through L , and (iv) es-caping through L . At the same time, we shall monitor theescape time (we will also use the terms escape rate and es-cape period) of the orbits.In Paper II we performed a similar orbit classificationhowever for energy levels relatively close to the escape en-ergy. In this case, we will proceed to energy levels a lothigher than the energy of escape and specifically in the in-terval E ∈ ( E L , − E > − x, z ) plane for four valuesof the energy is presented in Fig. 2(a-d). In these colour-coded grids we assign to each pixel a specific colour accord-ing to the corresponding type of the orbit. In particular,blue colour corresponds to regular non-escaping orbits, yel-low colour corresponds to trapped chaotic orbits, red colourcorresponds to orbits escaping through exit channel 1, whilethe initial conditions of orbits that escape through exit chan-nel 2 are marked with green colour. The outermost blacksolid line denotes the zero velocity curve which is defined as f ( x, z ) = Φ eff ( x, y = 0 , z ) = E. (7)For E = − .
26, which is an energy level just above theenergy of escape, we observe in Fig. 2a that more than halfof the ( x, z ) plane is covered by initial conditions of trappedchaotic orbits. The phenomenon of trapped chaos for energylevels above yet very close to the energy of escape has alsobeen observed in Paper II. Inside the vast chaotic sea wecan identify initial conditions of escaping orbits. Howeverescape is rather difficult for such energy levels. A substan-tial amount of the ( x, y ) plane is occupied by numerous is-lands corresponding to non-escaping regular motion. Thesestability islands correspond mainly to 1:1:0 resonant looporbits (see also Fukushige & Heggie 2000). As the value ofthe energy increases it is seen, in panels (b-d) of Fig. 2, thatthe portion of trapped chaotic orbits decreases, while at thesame time, the amount of escaping orbits grows rapidly. Onthe other hand, the structures of the stability islands re- mains almost unperturbed. In panel (d), where E = − . x, z ) plane is cov-ered by a high fractal mixture of escaping orbits. In thisdomain a high dependence of the escape mechanism on theparticular initial conditions of the orbits is observed. Thismeans that a minor change in the ( x , z ) initial conditionsof the orbit has as a result the star to escape through the op-posite escape channel which is of course a classical indicationof chaotic motion. It is interesting to note that especially atthe outer parts of the ( x, z ) plane, near the Lagrange points,basins of escape are present.The corresponding distribution of the escape time t esc of the orbits presented in Fig. 2(a-d) is illustrated in Fig.3(a-d). Once more, the scale of the escape time is revealedthrough a colour code. Specifically, initial conditions of fastescaping orbits are indicating with light reddish colors, whiledark blue/purple colors suggest high escape periods. Fur-thermore, initial conditions of non-escaping and trappedchaotic obits are shown in white. Note that the colour barhas a logarithmic scale. We observe in panels (a) and (b),where E = − .
26 and E = − .
24, respectively that the es-cape time of the orbits are huge corresponding to severalthousands of time units. However, as the value of the energyincreases the escape rates of the orbits are constantly beingreduced. This is true, because when E = − .
10 it is seenin panel (d) of Fig. 3 that the majority of the escaping or-bits need about a couple of hundred time units in order toescape. Looking carefully the time distributions it becomesevident that orbits with initial conditions in the fractal areasof the plot require a significant amount of time before theyescape from the cluster. On the contrary, the shortest escaperates have been measured for orbits inside the basins of es-cape where there is no dependence on the initial conditionswhatsoever.The escape dynamics of the dynamical system for fouradditional values of the energy is given in Fig. 4(a-d). Againdifferent colours are used for distinguishing between the fourtypes of the orbits (non-escaping regular, trapped chaotic,escaping through L and escaping through L ). In this caseall four energy levels are much higher than the energy ofescape. In Fig. 4a, where E = − .
5, we observe that escap-ing orbits cover about 80% of the ( x, z ) plane by formingwell-defined basins of escape. The size of stability islandscorresponding to non-escaping regular motion has been re-duced and they are mainly confined to the ( x <
0) part ofthe ( x, z ) plane. Therefore the majority of the ordered orbitsare in fact retrograde 1:1:0 resonant orbits (i.e., when a starrevolves around the star cluster in the opposite sense with re- We would like to emphasize that when it is stated that anarea is fractal we simply mean that it has a fractal-like geometrywithout conducting any specific calculations as in Aguirre et al.(2009). A local set of initial conditions that corresponds to a certainescape channel defines a basin of escape. The n : m : l notation we use for classifying the regular 3Dorbits is according to Carpintero & Aguilar (1998) and Zotos& Carpintero (2013). The ratio of those integers corresponds tothe ratio of the main frequencies of the orbit, where the mainfrequency is the frequency of the greatest amplitude for each co-MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 2.
Orbital structure of the ( x, z ) plane when (a-upper left): E = − .
26; (b-upper right): E = − .
24; (c-lower left): E = − . E = − .
10. The colour code is as follows: non-escaping regular orbits (blue), trapped chaotic orbits (yellow), escapingorbits through L (red), orbits escaping orbits through L (green). (For the interpretation of references to colour in this figure captionand the corresponding text, the reader is referred to the electronic version of the article.) spect to the motion of the cluster around the parent galaxy).Nevertheless, prograde 1:1:0 resonant orbits with x > E = − . ordinate. Main amplitudes, when having a rational ratio, definethe resonance of an orbit. E = − .
5, one may see that the stabilityislands are hardly visible, while escape regions dominate the( x, z ) plane. In particular, basins of escape corresponding to
MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters
Figure 3.
Distribution of the corresponding escape time t esc of the orbits on the ( x, z ) plane for the values of the energy presented inpanels (a-d) of Fig. 2, respectively. The darker the colour, the higher the escape time. Initial conditions of non-escaping regular orbitsand trapped chaotic orbits are shown in white. (For the interpretation of references to colour in this figure caption and the correspondingtext, the reader is referred to the electronic version of the article.) escape through L cover about 60% of the grid. Our compu-tations suggest that for E > − . x, z ) plane is E = − .
0. Our results are presented in panel (d) of Fig. 4. Forthis energy level we observe a very interesting phenomenon.Stability islands of non-escaping regular orbits emerge insidethe vast escape region in the left side ( x <
0) of the plane.Taking into account the outcomes shown in panels (a-c) ofFig. 4 it becomes more than evident that the area of thestability islands is constantly being reduced with increasing energy. However for E = − E = − E < −
1. In all previous cases, the vast majorityof non-escaping regular orbits corresponds to 1:1:0 resonantloop orbits, where the loop is parallel to the ( x, y ) plane (seepanel (a) of Fig. 5). On the other hand, when E = − MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 4.
Orbital structure of the ( x, z ) plane when (a-upper left): E = − .
5; (b-upper right): E = − .
0; (c-lower left): E = − . E = − .
0. The colour code is the same as in Fig. 2. (For the interpretation of references to colour in this figure captionand the corresponding text, the reader is referred to the electronic version of the article.)
The distribution of the escape time of the orbits on the( x, z ) plane for the values of the energy presented in Fig.4(a-d) is shown in Fig. 6(a-d), respectively. One may ob-serve that the results are very similar to those presentedearlier in Fig. 3(a-d), where we found that orbits with ini-tial conditions inside the basins of escape have the smallestescape rates, while on the other hand, the longest escapeperiods correspond to orbits with initial conditions either in the fractal regions of the plots, or near the boundaries of thestability islands of non-escaping regular motion.For the numerical integration of the initial conditionsof the three-dimensional orbits in each colour-coded grid onthe ( x, z ) plane, shown in Figs. 2 and 4, we needed about 1day of CPU time on a Quad-Core i7 2.4 GHz PC.Undoubtedly, it would be very illuminating to monitorthe evolution of the percentages of all types of orbits as a
MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters
Figure 5. (a-left): A 3D 1:1:0 loop orbit parallel to the ( x, y ) plane with initial conditions x = − . z = 0 . E = − . x = − . z = 0 .
25, for E = − .
0. In both cases y = p x = p z = 0,while the initial value of p y is obtained for the energy integral (5). function of the total orbital energy E . The correspondingdiagram is given in Fig. 7. We would like to point out thatfor constructing this diagram we used data from additionalcolour-coded grids apart from those given in Fig. 2 and Fig.4. We observe that just above the energy of escape trappedchaotic motion dominates as the corresponding initial con-ditions of the orbits cover about more than half of the ( x, z )plane. However their percentage drops very quickly as thevalue of the energy increases. In particular, for E > − E = − .
26. With increasing en-ergy the rate of non-escaping regular orbits decreases up to E = − .
5. At this energy level we observed the lowest rateof regular orbits which is about 1%. We could say that at E = − . x, y ) plane. For higher values of the energy ( E > − . x, y ) plane). The evolution of the escapingorbits is completely different with respect to the boundedorbits (regular or chaotic). Looking at Fig. 7 one can seethat for values of energy very close to the escape energy therates of both exit channels are equal which of course meansthat both channels are equiprobable. On the other hand, for E > − .
22 the rates of escaping orbits start to diverge. Be-ing more specific, the portion of escaping orbits through L increases more rapidly with respect to the portion of escap-ing orbit through L . Furthermore, for E > . L seems to saturates around 32%, whilethe rate of escapers through L continues to grow. At thehighest energy level studied, that is E = −
1, escaping orbitthrough exit channel 2 occupies about 63% of the ( x, z ) planewhich is about twice the rate of escapers through exit chan-nel 1. Therefore taking into account all the above-mentionedanalysis we may conclude that at low energy levels, wherethe fractality of the ( x, z ) plane is maximum, stars do notshow any particular preference regarding the escape chan-nels. On the contrary, at high enough energy levels, wherebasins of escape dominate, it seems that exit channel 2 (es-cape towards infinity) is twice more preferable with respectto exit channel 1 (escape towards the galactic center).Of course, the relative fraction of orbits escapingthrough saddle L and through saddle L strongly dependson our particular choice of the collection of initial conditions.The inversion symmetry of the whole system mentioned be-fore implies that to each orbit escaping through saddle L there is another orbit with inverted initial conditions whichescapes through saddle L . Therefore an integration overthe whole 5-dimensional energy shell of possible initial con-ditions would give equal escape rates through both saddles.And also a random collection of initial conditions from theinterior part of the 5-dimensional energy shell would leadto equal escape rates. For grids of initial conditions on anyparticular lower dimensional surface these escape rates de-pend on how this particular surface intersects the basins ofescape belonging to the two saddles.Before closing this section, we would like to present in MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 6.
Distribution of the corresponding escape time t esc of the orbits on the ( x, z ) plane for the values of the energy presented inpanels (a-d) of Fig. 4, respectively. The colour code is the same as in Fig. 3. (For the interpretation of references to colour in this figurecaption and the corresponding text, the reader is referred to the electronic version of the article.) Fig. 8 the evolution of the average value of the escape time < t esc > of the orbits as a function of the total orbital en-ergy E . We observe that very close to the critical energy ofescape the average escape time of the orbits correspond toabout 2000 time units. However, with increasing energy, therequired time for escape smoothly reduced until E = − < t esc > (cid:39) t esc continues but with a different slope. For E = − E = − .
26. Indeed, for E = − When we want to analyse and understand the dynamics of aHamiltonian system then in many cases it is easier to repre-sent the system by its Poincar´e map for fixed energy instead
MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters æææ æ æ æ æ æ æ æ æ æ æààà à à à à à à à à à àììì ì ì ì ì ì ì ì ì ì ìòòò ò ò ò ò ò ò ò ò ò ò - - - - - P e r c en t age H % L æ Trapped Chaotic à Non - Escaping Regular ì Escaping through L ò Escaping through L Figure 7.
Evolution of the percentages of all types of orbits onthe ( x, z ) plane, as a function of the total orbital energy E . Thedistribution shown applies to a particular plane of initial condi-tions. The asymmetry between the saddles L and L is causedby the choice of one branch of p y (i.e. one orientation of theinitial intersection of the plane y = 0). (For the interpretation ofreferences to colour in this figure caption and the correspondingtext, the reader is referred to the electronic version of the article.) - - - - - l og H < t e sc > L Figure 8.
Evolution of the logarithm of the average escape timeof the orbits (log ( < t esc > )), as a function of the total orbitalenergy E . of its flow. For a 2-dof Hamiltonian system the Poincar´emap acts on a 2-dimensional domain, therefore it can bepresented by 2-dimensional graphics and gives an instructiveoverview of the dynamics for fixed energy. Even better, if wehave the map for various values of the energy or of any otherparameter, then we see the development scenario of the sys-tem in a graphical form. For Poincar´e maps see the books ondynamical system theory Jackson (1991) and Lichtenberg &Lieberman (1993), while pictorial explanations can be foundin Abraham & Shaw (1992). In the present article we dealwith a Hamiltonian 3-dof system and then the Poincar´e mapfor fixed energy acts on a 4-dimensional domain and can nolonger be presented by 2-dimensional graphics. Then theimportant question is: How to gain an understanding of theimportant properties of this map? In addition, there is amore specific feature of these maps which is essential for allkinds of escape processes and which we like to transfer fromthe 2-dof to the 3-dof case. If the effective potential has anindex-1 saddle, then the flow over this saddle is the essen-tial step in the escape process and the escape is directedby some invariant subset sitting over this saddle and by itsstable and unstable manifolds.In 2-dof systems the invariant subset usually is an un-stable periodic orbit which in the Poincar´e map appears ashyperbolic fixed point. And the stable and unstable mani-folds of this periodic orbit in the flow (which appears as ahyperbolic fixed point in the map) direct and channel theflow over the saddle. The important point is that the 2-dimensional stable and unstable manifolds of the orbit inthe 3-dimensional flow for fixed energy (or the correspond-ing 1-dimensional stable and unstable manifolds of the fixedpoint in the 2-dimensional map) are of codimension 1 andtherefore they divide the phase space into regions of dif-ferent behaviour and are able to confine and channel themotion. The invariant subsets like the periodic orbit in the3-dimensional flow or the fixed point in the domain of themap are of codimension 2.When we go over to a 3-dof system and look for subsetsof analogous properties, then the codimensions and the nor-mal hyperbolicity are the essential properties which we mustmaintain. When we have some invariant subset with stableand unstable manifolds of codimension 1, then these mani-folds again divide the phase space and are able to confine andchannel the general motion. Accordingly, in a 4-dimensionalmap we need a 2-dimensional invariant subset sitting overthe index-1 saddle of the effective potential and having sta-ble and unstable manifolds of dimension 3. In addition, suchinvariant subsets should be persistent under general pertur-bations of the system in order to enable us to follow thescenario of the dynamics under parameter changes, for ex-ample under a change of the energy. Subsets having all theseproperties are known to mathematics, they are the normallyhyperbolic invariant manifolds (NHIMs) of codimension 2.They act as a kind of most important elements of the skele-ton of the whole dynamics. Additional general informationon NHIMs can be obtained in Wiggins (1994). In the presentsection we will study the NHIM sitting over the Lagrangepoint L of our system and in the following section we willstudy its implications for the tidal tails of the star cluster. MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
First let us consider the orbits which stay forever over anindex-1 saddle like the Lagrange point L . We will describethe important properties in a descriptive form without goinginto detailed calculations. The considerations of this subsec-tion are quite general and apply for any index-1 saddles ofthe effective potential in a rotating frame, for example theycoincide with the corresponding properties of the NHIM overthe Lagrange points for a barred galaxy as described in Pa-per III. In particular, the equations given in the beginning ofsection 4 and in the subsections 4.1 and 4.2 of Paper III canbe taken over for the present case immediately. Thereforewe do not repeat these equations here and we only give adetailed description of important properties. In the presentarticle we shall use the same notation as in Paper III, as faras possible.When the total orbital energy (numerical value of theHamiltonian in the rotating frame of reference) is above butclose to the saddle energy, then the dynamics over the sad-dle is well approximated by a quadratic expansion of theHamiltonian in the phase space coordinates. In this approx-imation the equations of motion are linear and accordinglythere are three normal modes of motion in the neighbour-hood of the saddle. Since we have an index-1 saddle, oneof these normal modes is unstable the other ones are sta-ble. The first stable mode is a vertical oscillation. When allavailable energy goes into this particular motion then wehave the vertical Lyapunov orbit, called Γ v in the follow-ing. The second stable mode is motion along a horizontalellipse. When all available energy goes into this particularmotion then we have the horizontal Lyapunov orbit, calledΓ h in the following. The third mode is unstable horizontalmotion which falls down exponentially along the unstabledirection of the effective potential. The general orbit staysin the neighbourhood of the saddle point for infinite timein the past and in the future if and only if it is a super-position of the two stable modes and does not contain anycontribution of the unstable motion of the third mode.For fixed total orbital energy there are two free pa-rameters for the general orbit which stays over the saddle.First, we can distribute the available energy between thetwo stable modes. Second, we can choose any phase shiftbetween these two stable modes. Accordingly, there is a 2-dimensional continuum of orbits which remain in the neigh-bourhood of the saddle for ever. This continuum forms a3-dimensional surface in the 5-dimensional flow for fixed en-ergy or a 2-dimensional surface in the 4-dimensional domainof the map. These surfaces are of codimension 2 and theyare unstable (hyperbolic) in the directions normal to the sur-face, i.e. these surfaces are the NHIMs we have been lookingfor. For energy sufficiently close to the saddle energy thetopology of the NHIM surface in the flow is the one of the3-dimensional sphere S . The NHIM of the map is repre-sented easiest as a curved disc. However, by contracting theboundary of the disc to a single point we can also representit as a 2-dimensional sphere S whenever this representationis more convenient.What happens under perturbations, for example whenthe energy increases and nonlinearities in the dynamics be-come important? Here one extremely important property of NHIMs becomes essential, namely their persistence un-der perturbations. With nonlinearities also the dynamics re-stricted to the invariant NHIM surface can develop instabil-ities. But for small perturbations the tangential instabilitiesare smaller than the normal instability and then the NHIMsurface survives as invariant surface. It may be shifted a lit-tle or deformed smoothly but it conserves its topology andconserves its NHIM properties. In particular, this impliesalso the persistence of its stable and unstable manifolds.Only for large perturbations the normal hyperbolicity canbe lost and then the NHIM may change qualitatively, maychange its topology or may decay and be lost. The reader canfind more information regarding the persistence propertiesin Berger & Bounemoura (2013); Eldering (2013); Fenichel(1971); Wiggins (1994), and for the bifurcations of NHIMsunder perturbations in Allahem & Bartsch (2012); MacKay& Strub (2014); Mauguiere et al. (2013); Teramoto et al.(2011, 2015a,b). Next we will investigate the perturbationscenario of the NHIM in our system under an increase of theenergy. The most prominent periodic orbits inside of the NHIM arethe Lyapunov orbits and some further periodic orbits relatedto the Lyapunov orbits. In the present subsection we shallstudy the development scenario of these periodic orbits un-der an increase of the energy as preparation of the study ofthe development scenario of the NHIM.In our system the vertical Lyapunov orbit Γ v followsan extremely simple scenario. The orbit Γ v is born at thesaddle energy E L and is stable in the directions tangentialto the NHIM and unstable in the directions normal to theNHIM. The qualitative stability properties do not changeunder an increase of E , only the numerical values of theeigenvalues change slowly. Accordingly, the orbit Γ v doesnot suffer any bifurcations along the whole scenario. Whenwe use the words tangential and normal here and in thefollowing then it always refers to the NHIM surface.Also the horizontal Lyapunov orbit Γ h is born at thesaddle energy as tangentially stable and normally unstable.Near the energy E = − .
275 it suffers a pitchfork bifurcationwhere it splits off a pair of tilted loop orbits, called Γ t andΓ t in the following. The orbit Γ v itself becomes tangen-tially unstable in this bifurcation. However, its tangentialinstability remains small compared to its normal instabil-ity. Accordingly, orbit Γ h remains part of the NHIM. Near E = − . h collides with anotherperiodic orbit called Γ s in the following which is normallyelliptic and the two colliding orbits disappear in a saddle-centre bifurcation in normal direction.The orbit Γ s itself is split off near E = − .
283 ina pitchfork bifurcation from another periodic orbit Γ r de-scribed below. The orbit Γ s is born stable in all directionsand under increasing energy it changes its stability prop-erties several times while it runs through several bifurca-tions where only two of them are of interest for us. Near E = − .
095 in a pitchfork bifurcation the orbit Γ s splitsoff two further tilted loop orbits which we call orbits Γ u and Γ u in the following. And finally, near E = − . MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters
Figure 9. (a-left): A collection of the most important periodic orbits in the configuration ( x, y, z ) space, for the energy E = − . x, y ), ( x, z ) and ( y, z ). The colour code is as follows: orbit Γ v (green), orbit Γ r (magenta), orbit Γ s (blue), orbit Γ h (red), orbits Γ t and Γ t (orange). (For the interpretation of references to colourin this figure caption and the corresponding text, the reader is referred to the electronic version of the article.) - - - - - E Figure 10.
Evolution of the x initial conditions of the periodicorbits as a function of the orbital energy E . The colour codeis the same as in Fig. 9, while orbits Γ u and Γ u are shown inbrown. (For the interpretation of references to colour in this figurecaption and the corresponding text, the reader is referred to theelectronic version of the article.) The periodic orbit Γ r mentioned above comes out of thepotential hole around the origin and it encircles the originin a shape which is symmetric under the reflection x → − x .For small values of the energy it is stable in all directions.Near E = − .
283 it suffers a pitchfork bifurcation where itbecomes unstable in one plane and splits off two new pe-riodic orbits where one of them is the orbit Γ s mentionedabove. The orbit Γ s is no longer invariant under the reflec-tion x → − x , however under this reflection the orbit Γ s ismapped into the second periodic orbit split off from the or-bit Γ r in the pitchfork bifurcation. Under increasing energyorbit Γ s moves in direction of increasing values of x , i.e. itmoves toward the saddle point L . The symmetry counter-part moves toward point L .The tilted loop orbits Γ t and Γ t are born normally un-stable and tangentially stable. They are part of the NHIM.Near E = − . t and Γ t collide with theother tilted loop orbits Γ u and Γ u respectively and allthese tilted loop orbits disappear in two symmetry relatedsaddle-centre bifurcations.In the following plots we present this scenario graphi-cally. Fig. 9(a-b) shows the mentioned periodic orbits in theposition space for the energy level E = − .
16. Note that forthis energy the orbits Γ u and Γ u do not yet exist. Panel (a)of Fig. 9 is a perspective view in the 3-dimensional ( x, y, z )position space, while panel (b) gives the 3 projections intothe 3 different coordinate planes. Green colour represents theorbit Γ v , magenta colour represents the orbit Γ r , blue colourrepresents the orbit Γ s , red colour represents the orbit Γ h ,while the two tilted loop orbits, Γ t and Γ t , are both plot-ted in orange colour. Note how the red and the blue curvealready become similar, so we can already anticipate the col- MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung - - - - - (a) - - - - - - - - (b) - - - - - (c) - - - - - - - - (d) Figure 11.
Evolution of (a upper left): the normal and (c lower left): the tangential traces of the monodromy matrix of the familiesof periodic orbits as a function of the total orbital energy E . Magnifications of panels (a) and (c) are shown in panels (b) and (d),respectively. The colour code is the same as in Fig. 10. The horizontal black dashed lines correspond to the critical values Tr crit = ± lision between these two orbits and their destruction in thesaddle-centre bifurcation under a small further increase ofthe energy. Also note how close the tilted loop orbits cometo the orbit Γ h in the projection into the ( x, y ) plane. Thisillustrates how the tilted loop orbits are split off from Γ h in a pitchfork in z -direction. We also see very well how theorbit Γ t is obtained from the orbit Γ t by a reflection in theplane z = 0.In Fig. 10 we show as function of the energy the x co-ordinate of all these orbits at the moment of intersection ofthe plane y = 0 in negative orientation. The use of colourscoincides with the one of Fig. 9, while the additional orbitsΓ u and Γ u are represented by the brown curve. Observethat a pair of tilted loop orbits gives a single curve only be-cause of the above mentioned symmetry under the reflection z → − z .The stability properties of the periodic orbits as func-tion of the energy are presented in Fig. 11(a-d). For allorbits, with the exception of the orbits Γ u and Γ u , the4-dimensional monodromy matrices have a natural decom-position into two 2-dimensional blocks and then we plot thetraces of these 2 × u and Γ u run through a complexspiralling case in the energy interval E ∈ [ − . , − . × × r is always ratherclose to +2. This orbit runs close to the origin basically in MNRAS465
Evolution of (a upper left): the normal and (c lower left): the tangential traces of the monodromy matrix of the familiesof periodic orbits as a function of the total orbital energy E . Magnifications of panels (a) and (c) are shown in panels (b) and (d),respectively. The colour code is the same as in Fig. 10. The horizontal black dashed lines correspond to the critical values Tr crit = ± lision between these two orbits and their destruction in thesaddle-centre bifurcation under a small further increase ofthe energy. Also note how close the tilted loop orbits cometo the orbit Γ h in the projection into the ( x, y ) plane. Thisillustrates how the tilted loop orbits are split off from Γ h in a pitchfork in z -direction. We also see very well how theorbit Γ t is obtained from the orbit Γ t by a reflection in theplane z = 0.In Fig. 10 we show as function of the energy the x co-ordinate of all these orbits at the moment of intersection ofthe plane y = 0 in negative orientation. The use of colourscoincides with the one of Fig. 9, while the additional orbitsΓ u and Γ u are represented by the brown curve. Observethat a pair of tilted loop orbits gives a single curve only be-cause of the above mentioned symmetry under the reflection z → − z .The stability properties of the periodic orbits as func-tion of the energy are presented in Fig. 11(a-d). For allorbits, with the exception of the orbits Γ u and Γ u , the4-dimensional monodromy matrices have a natural decom-position into two 2-dimensional blocks and then we plot thetraces of these 2 × u and Γ u run through a complexspiralling case in the energy interval E ∈ [ − . , − . × × r is always ratherclose to +2. This orbit runs close to the origin basically in MNRAS465 , 525–546 (2017) scape dynamics and NHIMs in star clusters
Figure 12.
Projections of the NHIM surfaces into the ( φ, L ) plane. The outermost solid closed curve corresponds to the horizontalLyapunov orbit. (a-upper left): E = −
3; (b-upper right): E = −
2; (c-lower left): E = − .
15; (d-lower right): E = − the almost rotationally symmetric potential of the interiorof the star cluster. This approximate rotational symmetryimplies approximate neutral stability under a rotation ofthe plane of this orbit and correspondingly it implies anapproximate value 2 of the trace. Also note that the stabilityproperties of a pair of tilted loop orbits coincide because oftheir interchange symmetry under the reflection z → − z . The NHIM in the map is a 2-dimensional invariant subsetof the domain of this map. This implies that the iteratedimages and preimages of any initial point from this subsetbelong to the same subset. Therefore the restriction of thePoincar´e map to this NHIM surface makes sense and it isa 2-dimensional Poincar´e map which can be represented by2-dimensional graphics. For previous examples of the con-struction and use of this restricted map see Gonzalez et al.(2014); Gonzalez & Jung (2015) and Paper III. Since NHIMsare persistent under perturbations also the existence of thismap is persistent under perturbations and is an interestingtool to study the perturbation scenario in 3-dof systems by2-dimensional graphics. It represents graphically the devel-opment scenario of the NHIM and it embeds the develop- ment scenario of the Lyapunov orbits studied in the previ-ous subsection into the development scenario of the wholeNHIM. For the idea how such a restricted map is constructednumerically see Gonzalez et al. (2014).For our numerical examples of Poincar´e maps we alwaysuse the intersection condition z = 0 with positive orienta-tion. As usual for the graphical presentation of Poincar´emaps we choose a moderate number of initial points in thedomain of the map and plot a large number of iterates ofthese initial points. The NHIM surface is a curved surfaceembedded into a 4-dimensional space. So to plot the mapon the surface we must use some kind of projection. In anal-ogy to Paper III it was instructive to use on one hand theprojection on the ( φ, L ) plane where φ = arctan( y/x ) and L = xp y − yp x and on the other hand perspective viewsof projections into the 3-dimensional ( x, y, p y ) space. Thisperspective view gives a suggestion of the curvature of theNHIM surface in the higher dimensional embedding space.For energy close to the saddle energy the dynamicscomes close to a linear dynamics and therefore the restrictedmap is similar to the Poincar´e map of a 2-dof anisotropicharmonic oscillator. As a numerical example of this case seepanel (a) of Fig. 12 which is constructed for E = − φ, L ) plane. The corresponding MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 13.
Projections of the NHIM surfaces into the ( x, y, p y ) phase space. The outermost solid closed curve corresponds to thehorizontal Lyapunov orbit. (a-upper left): E = −
3; (b-upper right): E = −
2; (c-lower left): E = − .
15; (d-lower right): E = − perspective plot in the ( x, y, p y ) space is shown in panel (a)of Fig. 13. The central fixed point on the NHIM representsthe Lyapunov orbit Γ v at the moment of its intersectionwith the plane z = 0 and the boundary of the domain is theLyapunov orbit Γ h which lies completely in the plane z = 0.If we want to represent also the orbit Γ h by a fixed pointof this map then we can contract the boundary to a sin-gle point. Thereby the domain of the map becomes a closedsurface with the topology of S . The rest of the domain ofthe map is filled by invariant lines which encircle the centralfixed point. They represent the quasi-periodic superpositionof vertical motion along the Lyapunov orbit Γ v and hori-zontal motion along the Lyapunov orbit Γ h . From the inside to the outside the energy in the vertical motion decreasesand the energy in the horizontal motion increases. For thecentral fixed point all energy is in the vertical motion andfor the boundary all energy is in the horizontal motion.When we increase the energy to the value E = − − , − E = − .
15. Now we arein the energy range of intermediate perturbation. In par-ticular near the boundary many secondary structures havebecome visible. In particular note the new periodic pointsnear φ = ± . L = − .
35. These two new structures be-
MNRAS465
MNRAS465 , 525–546 (2017) scape dynamics and NHIMs in star clusters long to the two tilted loop orbits Γ t and Γ t which are splitoff from the Lyapunov orbit Γ h , i.e. from the boundary, at E = − .
275 in a pitchfork bifurcation. Finally, in the panels(d) of Figs. 12 and 13 the energy has advanced to the value E = −
1. Now the Lyapunov orbit Γ h and the tilted looporbits have completely disappeared and the NHIM surfacehas lost its outer boundary.On the other hand, the inner part of the surface aroundthe fixed point belonging to the Lyapunov orbit Γ v remainsunaffected during the whole scenario. This is completely con-sistent with the property of the orbit Γ h not to suffer anybifurcations and to remain normally hyperbolic. This re-mains so up to high positive energies, and therefore also apart of the NHIM around the orbit Γ h remains for thesehigh energies.After having seen the development scenario of the Lya-punov orbits and of the NHIM in our present system, the fol-lowing questions naturally arise: Is this scenario well known?Does it fit into a more general scheme? How does it compareto the scenarios in similar systems?Unfortunately, the investigation of the general develop-ment scenarios of NHIMs is a problem little explored up tonow. The reader can find additional information on this in-teresting problem in the series of papers Allahem & Bartsch(2012); MacKay & Strub (2014); Mauguiere et al. (2013);Li et al. (2006); Teramoto et al. (2011, 2015a,b). Therefore,it is impossible to say, at this moment, whether or not ourpresent scenario fits into any general scheme. More specifi-cally we can compare our present scenario to the scenariosof other systems in a rotating frame of reference and hav-ing Lagrange points which are index-1 saddles of the effec-tive potential. Here 2 examples come to mind: First, therestricted three-body system as treated in Jorba & Masde-mont (1999), and second, the escape from a barred galaxyas it was described in Paper III. In Jorba & Masdemont(1999) the authors constructed the dynamics restricted tothe centre manifold of the saddle and studied the Poincar´emap of this restricted dynamics. Also they used the inter-section condition z = 0 for their map. The restriction ofthe dynamics to the centre manifold and the restriction ofthe dynamics to a NHIM are the same basic idea. Thus, wecan directly compare their restricted map to our restrictedmap. For instance, see Figs. 3 and 6 in Jorba & Masdemont(1999). Also in this system the vertical Lyapunov orbit doesnot suffer any bifurcation and the horizontal Lyapunov or-bit splits off a pair of tilted loop orbits. In this sense, thesystem of Jorba and Masdemont and the star cluster modelrun through the same basic scenario.In contrast, the scenario of the barred galaxy is qualita-tively very different. In the barred galaxy model the verticalLyapunov orbit runs through various bifurcations and splitsoff two different pairs of tilted loop orbits. In this systemthe horizontal Lyapunov orbit does not split off tilted looporbits. However, it splits off further horizontal orbits in apitchfork. Certainly the star cluster and the barred galaxybelong to clearly distinct scenarios. The best way to makeprogress in the extremely difficult but also important prob-lem of the development scenarios of NHIMs might be tostudy many examples, to recognize repeating patterns andto find phenomenological classifications. We hope that ourpresent example gives a useful contribution to this program. As mentioned above, the NHIMs have stable and unstablemanifolds of codimension 1 and these invariant manifolds ofthe NHIMs are responsible for the importance of the NHIMswith respect to the global dynamics of the system. If S E denotes the NHIM surface for the energy value E then we usethe symbols W s ( S E ) and W u ( S E ) for its stable and unstablemanifolds, respectively. W s ( S E ) and W u ( S E ) have each twobranches, one going outwards from the saddle and one goinginwards from the saddle respectively. W s ( S E ) contains suchorbits which converge towards S E in the future and W u ( S E )contains such orbits which converge towards S E in the past.This property gives already the idea how one can trace out W s ( S E ) and W u ( S E ) numerically. One simply takes a lot ofinitial conditions close to S E and let the orbits run. In thepast direction they approach automatically W s ( S E ) and inthe future they approach W u ( S E ).This is how panel (a) of Fig. 14 has been constructed forthe energy E = − . W s ( S E ) is plotted in green colour,while W u ( S E ) is plotted red. Of course, the manifolds livein the 5-dimensional energy shell of the phase space. Thefigure shows a perspective view of their projection into the3-dimensional position space. The energetically allowed re-gion of the position space in the inner part of the potentialis marked by gray colour. The figure shows the local partof the manifolds only, i.e. the orbits have been integratedfor a moderate time only. Globally, i.e. for large times, thesemanifolds grow complicated folds and tendrils. In panel (b)of Fig. 14 the projection of the stable and the unstable man-ifolds on the ( x, y ) plane are given. It is seen that the hori-zontal Lyapunov periodic orbit Γ h exist in the intersectionof the two manifolds. This periodic orbit is included into thesame figure as the blue solid curve.Fig. 14 shows the manifolds of the NHIM over the sad-dle point L only. Because of the inversion symmetry of thesystem mentioned before, there is an equal NHIM over thesaddle L . And one of these two NHIMs is obtained fromthe other one by a rotation around the z -axis by an angle π .Of course, the stable and unstable manifolds of one of theseNHIMs are also obtained from the stable and unstable man-ifolds of the other one by this symmetry operation. Globallythe stable manifolds of one NHIM intersect the unstablemanifold of the same NHIM and also the unstable manifoldof the other NHIM. Such intersections represent homoclinicand heteroclinic orbits which converge in the past as well asin the future to one of the NHIMs. Close to such intersec-tions exists also an infinity of periodic orbits which oscillatebetween the saddle points. General orbits move irregularlybetween these tangles.For an understanding of the escape process the follow-ing picture of the behaviour of typical orbits starting in theinner part of the potential is helpful. The orbit moves in theinner potential hole and at some time it finds itself closeto the local segment of the stable manifold of some sad-dle NHIM. Then it moves along this stable manifold to theneighbourhood of the NHIM and thereby also close to thesaddle. When the orbit does not start exactly on the stablemanifold then it stays in the neighbourhood of the NHIMfor a finite time only and leaves it again close to its unstablemanifold. It depends on fine details of the initial conditions MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 14. (a-left): The stable manifold W s ( S E ) (green) and the unstable manifold W u ( S E ) (red), when E = − .
24. The correspondingenergetically allowed region of the position space in the inner part of the potential is shown in transparent gray colour. (b-right): Theprojection of the stable and the unstable manifold on the configuration ( x, y ) space. The corresponding horizontal Lyapunov orbit Γ h isshown in blue. (For the interpretation of references to colour in this figure caption and the corresponding text, the reader is referred tothe electronic version of the article.) whether it leaves along the inner or the outer branch, i.e.whether it returns to the inner potential hole or escapes tothe outer region. The boundary between escape and returnis the stable manifold of the NHIM itself. By the mechanismjust described the NHIM is responsible for escape propertiesand for structures formed by the escaping orbits. More onthis point will be presented in the next section.In the previous subsection we studied the dynamics onthe NHIM surface itself. Now one can ask whether this re-stricted dynamics has any implications for the dynamics out-side of the NHIM. Here the stable and the unstable mani-folds of the NHIM are responsible for an affirmative answer.We can imagine W s ( S E ) and W u ( S E ) as a kind of cylindersover the NHIM which have an internal foliation accordingto the structures on the NHIM itself. And this foliation car-ries the structures found on the NHIM to far away regionsreached by W s ( S E ) and W s ( S E ). In this sense, the scenariofound in the restricted dynamics on the NHIMs is relatedto the dynamics in far away regions influenced by W s ( S E )and W u ( S E ). This section is devoted to the fate of escaping stars. In par-ticular, we will discuss what happens to stars when they passeither through L or L on their way to escape from the starcluster. It is well known that the majority of the stars of acluster escape from it, through the Lagrange points, with relatively low velocities . Escaping stars usually form ex-tended tidal tails (also known as tidal arms) (e.g., CapuzzoDolcetta et al. 2005; di Mateo et al. 2005; Just et al. 2009;K¨upper et al. 2008, 2010). Stars in tidal tails move under thegravitational influence of the axially symmetric potential ofthe parent galaxy due to the fact that the attraction of thestar cluster is practically negligible outside the tidal radius.Thus, due to the existence of Coriolis and centrifugal forcesescaping stars display a simple epicyclic motion along thetwo tidal tails.If there are stars in the interior region of the cluster withan energy high above the threshold energy E ( L ) then suchstars will leave the interior region fast and the interior regionwill have lost such stars long time ago. Let us now considerstars with an energy below the threshold but close to it andmoving in the inner part of the cluster. Such stars have oc-casional interactions among themselves and with other ob-jects and thereby their energy can be changed slightly andit may come a little above the threshold. Then such starsare exactly the ones for which the structure of the NHIMsover the Lagrange points L and L and their stable andunstable manifolds become highly relevant.First, these stars can come close to the saddle pointsalong the stable manifolds of the NHIMs, and then they Most stars escape as a result of two-body encounters. Theseescapers usually have an energy only slightly above the energy ofescape (Hayli 1970). Therefore, such escapers pass through oneof the two Lagrange points at slow speed.MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters - -
50 0 50y - -
505 10z
Figure 15.
The projection of the local segments of W u ( S E ) intothe ( x, y, z ) position space, when E = −
3. The manifolds fromthe NHIM over the saddle point L are shown in red colour, theones from the symmetrically placed NHIM over L are shown ingreen colour, while stars inside the interior region are shown inblue. (For the interpretation of references to colour in this figurecaption and the corresponding text, the reader is referred to theelectronic version of the article.) have two possibilities for their further motion. First, fromthe neighbourhood of the saddle they can return to the innerregion of the cluster along the inner branches of the unstablemanifolds of the NHIMs. Second, they can leave to the out-side of the cluster along the outer branches of the unstablemanifolds of the NHIMs. Which one of these two possibilitiesis realised depends on which side of the local segment of itsstable manifold the orbit approaches the NHIM. Stars whichreturn to the inside will come back to the neighbourhood ofa saddle point later and can eventually escape.When an orbit starts in the neighbourhood of the saddlepoint L then in forward direction (i.e. in the future) itconverges automatically against W u ( S E ). This observationprovides the idea for a numerical construction of W u ( S E ).We randomly select 1000 initial conditions close to L allwith the same energy E = − z =0. Approximately half of these orbits leave to the outsideimmediately along the outer branch of W u ( S E ). The otherhalf goes to the inside along the inner branch of W u ( S E ).Later after several revolutions in the inside this half of theorbits escapes partly over L and partly over L . In total allthese orbits together trace out the outer unstable manifoldsof both the NHIMs over L and L .Before looking at numerical plots we have to considerthree minor problems. First, it is too difficult to produce plots giving a good impression of these manifolds in the 5-dimensional energy shell of the phase space. Second and inaddition, our further discussions will focus on the positionspace. These points are taken care of by projecting W u ( S E )into the position space. Third, the complete unstable man-ifold has an infinite extension and has an infinity of foldsand tendrils and the outer branches leave to regions very farfrom the cluster. For our further discussion only the localsegments are important, i.e. the parts emanating directlyfrom the saddle region. In addition, we can expect that thetidal approximation no longer makes sense in regions farfrom the cluster. Therefore we cut off the unstable mani-fold, we restrict it by following the above mentioned orbitsover a finite time only. With these considerations in mind weplotted in Fig. 15 the projection into the ( x, y, z ) positionspace of the local segments of the unstable manifolds of thesaddle NHIMs. The time interval used for the cut off is [0,4.5]. In this figure, stars in the inner region of the clusterare coloured blue. Orbits leaving the cluster over the saddle L are plotted red, while orbits leaving over the saddle L are plotted green.This figure shows a very important property: The un-stable manifolds of the saddle NHIMs and therefore also thetidal tails of the cluster stay in a rather thin layer aroundthe plane z = 0. This means, that the tidal arms are alignedin the symmetry plane of the parent galaxy.Within the tidal approximation the two unstable man-ifolds, i.e. the one from saddle L and the one from saddle L , have the same shapes and are transformed into eachother by the inversion of horizontal coordinates. Withoutthe tidal approximation, i.e. under the full potential of theparent galaxy, these two manifolds should develop differentshapes on a large scale because they are developed underdifferent influences from the parent galaxy.For a better understanding of the behaviour of escap-ing stars, we reconfigured our numerical integration rou-tine for the computation of the basins of escape to yieldoutput of all orbits. For the initial condition of the orbitswe define a dense uniform three dimensional grid of size N x × N y × N z = 100 × × p x = p z = 0, whileboth branches of p y (obtained through the Jacobi integralof motion) are allowed. All the initial conditions of the or-bits lie inside the Lagrange radius of the star cluster, with E = −
3. We numerically integrate the initial conditions ofthe three-dimensional orbits and we record the output of allorbits. This allow us to monitor the formation as well asthe time-evolution of the tidal tails constructed by the starsthat escape through L and L .Initially, at t = 0 all orbits are regularly distributedwithin the Lagrange radius. In Fig. 16(a-d) we present thetime-evolution of the ( x, y, z ) position of the stars. It is seenthat as time evolves both tidal tails (leading and trailing)grow along the y -axis however, the values of the x and z co-ordinates remain always very low ( | x | < | z | < . t max = 30 so that y max (cid:39) r L . We observe in panel(a) of Fig. 16 that at t = 7 . L or L . As time goes by we observe in panels (b-d) thatthe two tidal tails grow in size and at t = 30 time units thestellar structures are fully developed, while at the same time MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
Figure 16.
The distribution of the position of stars in the configuration ( x, y, z ) space initiated ( t = 0) within the Lagrange radius, for E = −
3. The leading tail (red) contains stars that escaped through L , while the trailing tail (green) contains stars that escaped through L . Stars moving in bounded (trapped chaotic or non-escaping regular) motion inside the Lagrange radius are shown in blue. Note thatthe development of the tidal tails takes place very close to the primary ( x, y ) plane. (For the interpretation of references to colour in thisfigure caption and the corresponding text, the reader is referred to the electronic version of the article.) the area inside the tidal radius is uniformly filled. Here itshould be pointed out that the standard integration code,which was used to obtain the results presented in Figs. 2 and4, records a point at every step of the numerical integration.In Fig. 16(a-d) on the other hand, the density of points along one star orbit is taken to be proportional to the velocity ofthe star. Being more specific, a point is plotted (showingthe position of a star), if an integer counter variable whichis increased by one at every integration step, exceeds thevelocity of the star. Following this technique we can simu- MNRAS465
3. The leading tail (red) contains stars that escaped through L , while the trailing tail (green) contains stars that escaped through L . Stars moving in bounded (trapped chaotic or non-escaping regular) motion inside the Lagrange radius are shown in blue. Note thatthe development of the tidal tails takes place very close to the primary ( x, y ) plane. (For the interpretation of references to colour in thisfigure caption and the corresponding text, the reader is referred to the electronic version of the article.) the area inside the tidal radius is uniformly filled. Here itshould be pointed out that the standard integration code,which was used to obtain the results presented in Figs. 2 and4, records a point at every step of the numerical integration.In Fig. 16(a-d) on the other hand, the density of points along one star orbit is taken to be proportional to the velocity ofthe star. Being more specific, a point is plotted (showingthe position of a star), if an integer counter variable whichis increased by one at every integration step, exceeds thevelocity of the star. Following this technique we can simu- MNRAS465 , 525–546 (2017) scape dynamics and NHIMs in star clusters late, in a way, a real N -body simulation of the tidal tailsevolution of the star cluster, where the density of stars willbe highest where the corresponding velocity is lowest. Addi-tional numerical simulations, not presented here, reveal thatthe same behavior occurs for other (lower or higher) valuesof energy.With a much closer look at Fig. 16(a-d) it is easy forsomeone to observe that the two developed tidal tails arenot symmetrical. Our calculations indicate that during theentire time interval of the development of the stellar struc-tures the portion of stars that escape through L is slightlyelevated with respect to the amount of stars that escapethrough L . Being more precise, exit channel 2 appears tobe about 6% more preferable than exit channel 1. This ob-servation comes in an agreement with the results presentedearlier in Fig. 7, where we seen that the percentage of es-caping orbits through L is always higher than that the rateof escaping orbits through L . Only when the energy levelis relatively close to the critical energy of escape these twopercentages coincide due to the complete fractality of thephase space. At this point, we would like to clarify that thedeviation of the escape rates of the two channels is due tothe particular choice of the initial conditions of the orbits. The aim of this work was to numerically investigate the3-dof escape dynamics of a star cluster embedded in thesteady tidal field of a parent galaxy. We considered energylevels larger than the escape energy where the equipoten-tial surface opens and two exit channels appear, throughwhich the test particles (stars) can escape to infinity. Wemanaged to distinguish between ordered versus chaotic andtrapped versus escaping orbits. We also located the basinsof escape leading to different exit channels, finding at thesame time correlations with the corresponding escape timeof the orbits. Our extensive and thorough numerical explo-ration strongly suggests that the overall escape mechanismis a very complicated procedure.The NHIMs over the saddle points act as the most im-portant elements of the skeleton of the global dynamics asthey direct the escape process over the saddles. Therefore westudied the development scenario as function of the total or-bital energy for the NHIM over the saddle point L in detailand this includes the development scenario of the Lyapunovorbits contained in this NHIM. The main tool to study theNHIM has been the numerical construction of the restrictedPoincar´e map on the NHIM. It provides simultaneously theinvariant surface itself and the dynamics on it.The stable and the unstable manifolds of the NHIMsplay an essential role because they transport important dy-namical structures created on the NHIM to far away regionsof the phase space and they channel and direct the escapeflow over the saddles. As a consequence, the escaping orbitsfollow the unstable manifolds of the NHIMs and thereby theprojection into the position space of these unstable mani-folds delimits the tidal tails (or arms) of the cluster formedby the escaping stars.The main numerical results of our research can be sum-marized as follows:(i) Regions of bounded regular motion were found to co- exist with extended basins of escape composed of initialconditions of orbits that escape through the exit channelslocated near the two saddle points.(ii) Bounded regular orbits mainly correspond to simple1:1 loop three-dimensional orbits, while other secondary res-onant orbits were also observed. The shift of the value ofthe total orbital energy mostly influences the orientation ofthese loop orbits. In particular, we found that at relativelyhigh energy levels the 1:1 loop orbits become inclined to the( x, y ) plane.(iii) At energy levels above, yet very close to the criticalenergy of escape, we detected the existence of a substantialamount of chaotic orbits which remain trapped inside theLagrange radius for vast time intervals. However the per-centage of these trapped chaotic orbits heavily reduces aswe proceed to higher energy levels.(iv) A strong correlation between the value of the Jacobiintegral of motion and the extent of the basins of escapewas found to exist. More precisely, for low energy levels thestructure of the ( x, z ) plane exhibits a large degree of fractal-ization which implies that the majority of the orbits escapechoosing randomly escape channels. As the value of the en-ergy increases the orbital structure becomes less and lessfractal, while several basins of escape emerge and dominate.(v) It was revealed that the percentages of the orbits thatescape through the two escape channels are almost equalonly very close to the energy of escape, where the degreeof fractalization is still strong. For higher values of the en-ergy on the other hand, their rates start to diverge and atrelatively high energy levels channel 2 (around L ) seemsto be two times more preferable with respect to channel 1.However this behaviour is due to the particular choice of theinitial conditions of the orbits.(vi) Our numerical analysis on the development scenarioregarding the families of the periodic orbits contained in theNHIMs suggests that the orbital content of the network ofthe periodic orbits is much simpler with respect to otherphysical systems (i.e., a 3-dof model of a barred galaxy). Inother words, our computations indicate that there are lessbifurcated orbits of the main Lyapunov orbits.(vii) We provided numerical evidence that stars with ini-tial conditions in the close vicinity of the NHIMs of the twoLagrange points L and L can form tidal tails (or tidalarms) upon escape through the exit channels around thesaddles.We hope that the results of the present numerical ex-ploration shed some light on the properties of the escapedynamics as well as on the role of the normally hyperbolicinvariant manifolds (NHIMs) in tidally limited star clusters. ACKNOWLEDGMENTS
One of the authors (CJ) thanks DGAPA for financial sup-port under grant number IG-100616. We would like to ex-press our warmest thanks to the anonymous referee for thecareful reading of the manuscript and for all the apt sugges-tions and comments which allowed us to improve both thequality and the clarity of our paper.
MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung
REFERENCES
Abraham R.H., Marsden J.E., 1978, Foundations of Mechan-ics 2nd ed., The Benjamin Cummings Publishing Company,Reading MassachusettsAbraham R.H., Shaw C.D., 1992, Dynamics: The Geometry ofBehavior, Addison-Wesley, Redwood City, CA, 2nd ed.Aguirre J., Viana R.L., Sanju´an M.A.F., 2009, Rev. Mod. Phys.,81, 333Allahem A., Bartsch T., 2012, J. Chem. Phys., 137, 214310Ambartsumian V.A., 1938, Ann. Leningrad State Univ., 22, 19Baumgardt H., 2001, MNRAS, 325, 1323Belokurov V., Evans N.W., Irwin M.J., Hewett P.C., WilkinsonM.I., 2006, ApJ, 637, L29Berger P., Bounemoura A., 2013, Dynamical Systems, 28, 567Binney J., Tremaine S., 2008, Galactic Dynamics, Princeton Univ.Press, Princeton, USABressert E., Bastian N., Gutermuth R., Megeath S.T., Allen L.,et al., 2010, MNRAS, 409, L54Capuzzo Dolcetta R., Di Matteo P., Miocchi P., 2005, AJ, 129,1906Carpintero D.D., Aguilar L.A. 1998, MNRAS, 298, 1Chandrasekhar S., 1942, Principles of Stellar Dynamics. Univ.Chicago Press, ChicagoChoi J.-H., Weinberg M.D., Katz N., 2007, MNRAS, 381, 987Dehnen W., Odenkirchen M., Grebel E.K., Rix H.-W., 2004, AJ,127, 2753di Matteo P., Capuzzo Dolcetta R., Miocchi P., 2005, Cel. Mech.Dyn. Astron., 91, 59Eldering J., 2013, Normally Hyperbolic Invariant Manifolds, TheNoncompact Case, Atlantis PressErnst A., Just A., Spurzem R., Porth O., 2008, MNRAS, 383, 897(Paper I)Feast M., Whitelock P., 1997, MNRAS, 291, 683Fellhauer M., Evans N.W., Belokurov V., Wilkinson M.I., GilmoreG., 2007, MNRAS, 380, 749Fenichel N., 1971, Indiana Univ. Math. J., 21, 193Fukushige T., Heggie D.C., 2000, MNRAS, 318, 753Gonzalez F., Drotos G., Jung C., 2014 J. Phys. A: Math. Theor.,47, 045101Gonzalez F., Jung C., 2015, J. Phys. A: Math. Theor., 48, 435101Grillmair C.J., Freeman K.C., Irwin M., Quinn P.J., 1995, AJ,109, 2553Hayli A., 1970, A&A, 7, 17H´enon M., 1960, Ann. Astrophys., 23, 668Jackson E.A., 1991, Perspectives of nonlinear dynamics, Cam-bridge University PressJohnston K.V., Sigurdsson S., Hernquist L., 1999, MNRAS, 302,771Jorba A., Masdemont J., 1999, Physica D, 132, 189Jung Ch. Pott S., 1989, J. Phys. A22, 2925Jung Ch. Scholz H.-J., 1988, J. Phys. A 21, 2301Jung Ch., Taylor H.S., 1994, Z. Phys. D, 32, 79Jung Ch., Zotos E.E., 2015, PASA, 32, e042Jung Ch., Zotos E.E., 2016, MNRAS, submitted, (Paper III)Just A., Berczik P., Petrov M., Ernst A., 2009, MNRAS, 392, 969Kharchenko N., Scholz R.-D., Lehmann I., 1997, A&AS, 121, 439King I.R., 1959, AJ, 64, 351King I.R., 1962, AJ, 67, 471King I.R., 1966, AJ, 71, 64Koch A., Grebel E.K., Odenkirchen M., Mart´ınez-Delgado D.,Caldwell J.A.R., 2004, AJ, 128, 2274Kollmeier J.A., Gould A., 2007, ApJ, 664, 343Kruijssen J.M.D., 2012, MNRAS, 426, 3008K¨upper A.H.W., Macleod A., Heggie D.C., 2008, MNRAS, 387,1248K¨upper A.H.W., Kroupa P., Baumgardt H., Heggie D.C.,Lada C.J., Lada E.A., 2003, ARA&A, 41, 57 Lee K.H., Lee H.M., Fahlman G.G., Lee M.G., 2003, AJ, 126, 815Lee K.H., Lee H.M., Sung H., 2006, MNRAS, 367, 646Lehmann I., Scholz R.-D., 1997, A&A, 320, 776Li C.B., Shoujiguchi A., Toda M., Komatsuzaki T., 2006 Phys.Rev. Lett., 97, 028302Li Y., Luo A., Zhao G., Lu Y., Ren J., Zuo F., 2012, ApJ, 744,L24Lichtenberg A.J., Lieberman M.A., 1993, Regular and stochasticmotion, New York: Springer VerlagLyapunov A.M., 1907, Ann. Fac. Sci. Toulouse 9, 203MacKay R.S., Strub D.C., 2014, Nonlinearity, 27, 859Marsden J.E., 1992, in Lectures on mechanics, London Math. Soc.Lecture Note Series, Vol. 174, Cambridge University PressMauguiere F.A.L., Collins P., Ezra G.S., Wiggins S., 2013 Int. J.Bif. Chaos, 23, 1330043Montuori M., Capuzzo-Dolcetta R., di Matteo P., Lepinette A.,Miocchi P., 2007, ApJ, 659, 1212Odenkirchen M., Grebel, E.K., Rockosi, C.M., Dehnen, W., Ibata,R., et al., 2001, ApJ, 548, L165Press H.P., Teukolsky S.A, Vetterling W.T., Flannery B.P., 1992,Numerical Recipes in FORTRAN 77, 2nd Ed., CambridgeUniv. Press, Cambridge, USARockosi C.M., Odenkirchen, M., Grebel, E.K., Dehnen, W., Cud-worth, K.M., Gunn, J.E., et al., 2002, AJ, 124, 349Ross D.J., Mennim A., Heggie D.C., 1997, MNRAS, 284, 811Royer F., 1997, in R.M. Bonnet, E. Høg, P.L. Bernacca, L. Emil-iani, A. Blaauw, C. Turon, J. Kovalevsky, L. Lindegren, H.Hassan, M. Bouffard, B. Strim, D. Heger, M.A.C. Perryman& L. Woltjer ed., Hipparcos - Venice ’97 Vol. 402 of ESASpecial Publication, Populations among High-Velocity Early-Type Stars. pp 595-598Silva M.D.V., Napiwotzki R., 2011, MNRAS, 411, 2596Skokos C., 2001, Journal of Physics A, 34, 10029Smale S., 1970a, Inv. Math. 10, 305Smale S., 1970b, Inv. Math. 11, 45Spitzer L., 1940, MNRAS, 100, 396Teramoto H., Toda M., Komatsuzaki T., 2011, Phys. Rev. Lett.,106, 054101Teramoto H., Toda M., Komatsuzaki T., 2015a, Nonlinearity, 28,2677Teramoto H., Toda M., Takahashi M., Kono H., Komatsuzaki T.,2015b, Phys. Rev. Lett., 115, 093003Wiggins S., 1994, Normally Hyperbolic Invariant Manifolds inDynamical Systems, Berlin: Springer VerlagWolfram S., 2003, The Mathematica Book. Wolfram Media,ChampaignYim K.-J., Lee H.M., 2002, JKAS, 35, 75Zotos E.E., 2015a, MNRAS, 446, 770 (Paper II)Zotos E.E., 2015b, MNRAS, 452, 193Zotos E.E., Carpintero D.D., 2013, CeMDA, 116, 417
APPENDIX A: RELATION BETWEEN THETOTAL AND THE EFFECTIVE POTENTIAL
In previous works on the subject of tidally limited star clus-ters (e.g., Papers I and II) the effective potentialΦ eff ( x, y, z ) = Φ cl ( x, y, z ) + 12 (cid:0) κ − (cid:1) x + 12 ν z , (8)has been used. However in this paper we decided to use thetotal potential Φ t given in Eq. (2). This is necessary becausewe used the canonical form instead of the velocity form forthe equations of motion and the variational equations. Inthe following, we shall justify our choice and we will alsomake clear the connection between them. For any system ofcoordinates in position space we can write the equations of MNRAS , 525–546 (2017) scape dynamics and NHIMs in star clusters motion either in velocity form or in the canonical form. Thecanonical form always has the advantage that the symplecticstructure of Hamiltonian dynamics becomes evident. Thenwe can apply the whole mathematical machinery of symplec-tic differential geometry without unnecessary complications.For a mathematically oriented presentation of Hamiltoniandynamics see Abraham & Marsden (1978).Hamiltonian systems have some extra structure com-pared to more general dynamical systems, it is the sym-plectic structure. However, this symplectic structure onlybecomes evident in canonical coordinates. These are coor-dinates where the symplectic differential form ω has thesimple functional form according to the Darboux theorem ω = N (cid:80) j =1 dp j ∧ dq j , where j runs over all degrees of freedom.These coordinates are not unique. A transformation fromone coordinate system with the Darboux form to anotherone which also has the Darboux form is called canonicaltransformation. Only in these coordinates the area and vol-ume conservation in various dimensions of the flow and ofthe Poincar´e map is obvious. It is trivial to show that thetime evolution itself is a canonical transformation. And, im-portant for the following, only in canonical coordinates therelation between the total potential Φ t and the effective po-tential Φ eff is easy to understand.Let the Hamiltonian in the inertial frame have the stan-dard form H = T + Φ t , where T is the usual kinetic energyand Φ t is the total potential as a function of position co-ordinates q (actually q is a multi component). The angularmomentum is the generating function for rotations. Hereone important property of canonical coordinates comes in:Any quantity is the generating function for translations inthe conjugate direction. Therefore the transition to a ro-tating coordinate system with rotation frequency Ω is justthe coordinate system where φ moves with constant veloc-ity and therefore it is just generated by − Ω L , where L isthe angular momentum. Therefore in order to obtain theHamiltonian in the rotating frame of reference we just add − Ω L to the Hamiltonian in the inertial frame. The effectivepotential Φ eff is defined as Φ eff ( q ) = min p H ( q, p ), wherewe have to take the minimum of H over all possible val-ues of p for fixed q . When H is a smooth function thenthe minimum is obtained at the point where ∂H/∂p = 0.But this partial derivative is just the velocity according tothe Hamiltonian equations of motion. In our particular casethis condition of velocity 0 means that p x = − Ω y , p y = Ω x and p z = 0. Thus we obtain the effective potential when weplug in these values for the momenta into the Hamiltonianin the rotating frame of reference. In our case the result isΦ eff = Φ t − Ω (cid:0) x + y (cid:1) . More details regarding the generalmethod for constructing effective potentials can be found inJung & Taylor (1994); Marsden (1992); Smale (1970a,b).Now let us imagine that we have a good approximationfor Φ eff , usually it is some truncation of a power series ex-pansion. However in H we need Φ t . Therefore we invert theabove relation and we write Φ t = Φ eff + Ω (cid:0) x + y (cid:1) andfinally we insert the approximation for Φ eff into the Hamil-tonian. For the tidally limited cluster we are exactly in thesame situation. The tidal approximation gives a simple func-tional form for Φ eff and we can apply all this machinery tothe star cluster. APPENDIX B: LOOPS OF ORBITS AROUNDTHE NHIM
If an initial condition lies exactly on the stable manifold ofsome unstable (hyperbolic) invariant subset, then the corre-sponding orbit converges to this subset and remains on it forever. Of course, numerically we can never start exactly on astable manifold, we can only start close to it. Then the orbitfirst approaches the invariant subset, stays close to it for afinite time and later leaves the neighbourhood of this subsetagain near and along its unstable manifold. The time whichthe orbit spends in the neighbourhood of the invariant sub-set depends on how close to the stable manifold the orbithas been started. If we come closer to the stable manifoldby a factor eigenvalue of the invariant subset, then the orbitmakes one more loop close to the invariant subset.For 2-dof systems the relevant subsets are unstable (hy-perbolic) periodic orbits. And then coming closer to its sta-ble manifold by one eigenvalue means making one additionalrevolution close to this periodic orbit. For an instructivenumerical example of this behaviour see Fig. 7 in Jung &Scholz (1988). Additional and more general explanations ofthe scaling behaviour of escaping orbits in the neighbour-hood of stable manifolds of invariant subsets can be foundin Jung & Pott (1989).The above-mentioned early examples treated 2-dof sys-tem. However, the basic idea holds the same for any numberof degrees of freedom. For NHIMs of codimension 2 the rel-evant eigenvalue is the one in normal direction. A numericalexample from the present system of the star cluster is pre-sented in Fig. 17. From one part of the figure to the nextone 2 more decimal digits have been included into the ini-tial condition of the x coordinate, And the sequence of initialconditions converges to the stable manifold of the NHIM. Weclearly see how the number of loops which the orbit makesover the NHIM, and also around the Lagrange point L ,increases with increasing precision of this initial condition. This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS , 525–546 (2017) E. E. Zotos & Ch. Jung - - - - - L (a) - - - - - L (b) - - - - - L (c) - - - - - L (d) - - - - - L (e) - - - - - L (f ) Figure 17.
The projection on the ( x, y ) plane of a three dimensional orbit with initial conditions: y = 0 , z = 0 . , p x = p z = 0. Thevalue of p y is obtained from the energy integral with E = − .
5. The numerical value of the initial x coordinate is x = 0 . x varies through the panels as follows: (a): four digits; (b): six digits; (c): eight digits; (d):ten digits; (e): twelve digits; (f): fourteen digits. The position of the Lagrange point L is indicated by a blue dot, while the outermostred solid line is the horizontal Lyapunov orbit. We observe that the number of loops around L increases with increasing accuracy of x .(For the interpretation of references to colour in this figure caption and the corresponding text, the reader is referred to the electronicversion of the article.) MNRAS465