Practical application of KAM theory to galactic dynamics: II. Application to weakly chaotic orbits in barred galaxies
MMNRAS , 1–17 (2015) Preprint 17 September 2018 Compiled using MNRAS L A TEX style file v3.0
Practical application of KAM theory to galactic dynamics:II. Application to weakly chaotic orbits in barred galaxies
Martin D. Weinberg (cid:63)
Department of AstronomyUniversity of Massachusetts, Amherst MA 01003-9305, USA
17 September 2018
ABSTRACT
Owing to the pioneering work of Contopoulos, a strongly barred galaxy is known tohave irregular orbits in the vicinity of the bar. By definition, irregular orbits can notbe represented by action-angle tori everywhere in phase space. This thwarts pertur-bation theory and complicates our understanding of their role in galaxy structure andevolution. This paper provides a qualitative introduction to a new method based onKAM theory for investigating the morphology of regular and irregular orbits based ondirect computation of tori described in Weinberg (2015a) and applies it to a galaxydisc bar. Using this method, we find that much of the phase space inside of the bar ra-dius becomes chaotic for strong bars, excepting a small region in phase space betweenthe ILR and corotation resonances for orbits of moderate ellipticity. This helps explainthe preponderance of moderately eccentric bar-supporting orbits as the bar strengthincreases. This also suggests that bar strength may be limited by chaos! The chaos re-sults from stochastic layers that form around primary resonances owing to separatrixsplitting. Most investigations of orbit regularity are performed using numerical com-putation of Lyapunov exponents or related indices. We show that Lyapunov exponentspoorly diagnose the degree of stochasticity in this problem; the island structure in thestochastic sheaths allow orbit to change morphology while presenting anomalouslysmall Lyapunov exponent values (i.e. weak chaos). For example, a weakly chaoticorbit may appear to change its morphology spontaneously, while appearing regularexcept during the change itself. The numerical KAM approach sensitively detectsthese dynamics and provides a model Hamiltonian for further investigation. It mayunderpredict the number of broken tori for strong perturbations.
Key words: methods: numerical — stars: kinematics and dynamics — galaxies:kinematics and dynamics, formation, structure
The dynamics of barred galaxies has influenced our under-standing of secular evolution and the importance of clas-sical indeterminism or chaos over the last several decades(see Contopoulos 2002). In the standard non-linear dynam-ics literature, chaos is often classified by the magnitude andcharacter of the exponential divergence of nearby trajecto-ries in phase space. If the chaos is strong, one may measureexponential divergence in at least one dimension, and thismay result in orbital diffusion in otherwise conserved quan-tities of the motion. Conversely, some orbits may remainin a small region of phase confined by regions of regular-ity, possibly with small or undetectably positive exponen-tial divergence; this chaos is often called ‘weak’. Athanas- (cid:63)
E-mail: [email protected] soula, Romero-G´omez & Masdemont (2009); Athanassoulaet al. (2009); Manos & Athanassoula (2011); Contopoulos &Harsoula (2010); Kaufmann & Contopoulos (1996); Patsis,Athanassoula & Quillen (1997); Romero-G´omez et al. (2006,2007); Manos (2008); Tsoutsis et al. (2009); Brunetti, Chi-appini & Pfenniger (2011); Bountis, Manos & Antonopou-los (2012) have all discussed the role of chaos in barredgalaxies. In addition, (Cachucho, Cincotta & Ferraz-Mello2010; Giordano & Cincotta 2004; Manos & Machado 2014)have recently discussed to the implications of irregularityin n-body simulations of barred galaxies. Clearly, it is well-established that chaos is pervasive in the presence of strongpatterns. The work described here is motivated by the de-sire to connect the onset of chaos to methods of classicalperturbation theory that underpins our theory of the strongpatterns themselves. The term chaos has broad context andappears to have many definitions. For the remainder of this c (cid:13) a r X i v : . [ a s t r o - ph . GA ] A ug M. D. Weinberg paper, we will call any orbit that can not be fully describedby its conserved quantities as irregular or chaotic synony-mously, and we will refer to periods of exponential sensitivityto initial conditions as stochastic .The most commonly used technique to diagnose chaosin time-independent potentials is Lyapunov exponent anal-ysis. These exponents, λ measure the exponential diver-gence of two initially close initial conditions | X (0) − X (0) | = O ( (cid:15) ) in phase space in the form | X ( t ) − X ( t ) | =exp( λt ) | X (0) − X (0) | . Lyapunov exponent analysis is ro-bust: it can be applied to any equations of motion regardlessof complexity, although those that can be linearised ana-lytically have some computational advantage. Let us con-sider the application of this method to a background modeldescribed by its full complement of actions, I , and a per-turbation which induces irregularity. Chaotic orbits by na-ture appear to diffuse through the previously invariant toridescribed by I and may be detected by positive Lyapunovexponents. Other chaos detection schemes, such as Fourierspectral analysis, short-term Lyapunov analysis, and the re-cent the Generalised Alignment Index (GALI) may be usedin more general time-dependent systems, but we will notconsider these here. Although many of these have been pro-ductively and creatively used to investigate astronomicalsystems, these exponents and indices can be difficult to linkto the underlying dynamical mechanisms relevant for as-tronomical problems. For example, Lyapunov exponents forweak chaos may be small or zero, but induced change in theorbits may be qualitatively large and result in a completelydifferent orbit morphology. Weak chaos appears to dominatefor barred systems and we will illustrate the effect on specificorbits in section 4.In this paper, we will explore an alternative approachthat reconstructs the perturbed tori directly using a numer-ical technique based on Kolmogorov-Arnold-Moser (here-after, KAM) theory, described in Weinberg (2015a, Paper1). The technique attempts to find a new torus (regularorbit) using the constructive procedure outlined in manyproofs of the KAM theorem. This method has several dis-tinct advantages: 1) it is computationally efficient relativeto exhaustive numerical integration of orbits; 2) it providesa statistical assessment of the regular-orbit fraction in ac-tion space of the unperturbed problem; this in turn, allowsa direct connection with our standard dynamical descrip-tion of patterns and instabilities in galactic dynamics; 3) itprovides self-diagnostic information about the broken torus,e.g. which resonances are responsible. Our method does notfollow the KAM procedure precisely. Since we are interestedin locating and diagnosing broken tori, we do not attemptto find nearby action values that yield valid tori, but ratherinterpret the failure in a statistical sense. Although the pro-posed approach is perturbative, it applies to both generaltime-dependent and time-periodic perturbations. Paper 1suggests that the two tools together give a more completepicture than either one alone. We will see that this is true forthe present application as well. Our main goal is to charac-terise the use of this new method on a popular astrophysicaldynamics problem and compare with the more commonlyused Lyapunov exponent method.For completeness, we will begin with a short review ofboth Lyapunov exponent analysis and the numerical KAMapproach. A more in-depth but non-mathematical descrip- tion of the KAM-motivated ideas can be found in section 2.Consider a general n-dimensional dynamical system d X dt = F ( X , t ) , (1)where X is an n -dimensional vector (e.g. 6-dimensionalphase space) and F is smooth continuous vector function(e.g. a physical force function). Our initial condition is givenby X (0) at t = 0. The n Lyapunov exponents of the systemare defined by the rate of the logarithmic increase of theaxes of an infinitesimal sphere of states around X (0). Tocompute this rate, one may use the tangent map given bythe first-order expansion of equation (1): ∂δ X dt = ( ∇ · F ) · δ X (2)where ∇· F is the n × n Jacobian matrix ( ∇· F ) ij ≡ ∂F i /∂X j .One of the standard methods used to determine the full Lya-punov spectrum is due to Shimada & Nagashima (1979) andBenettin et al. (1980) who use a Gram-Schmidt reorthonor-malisation procedure to identify the principle axes of theexponentially diverging and converging sphere. An explicitsource code for computations based on this procedure isgiven by Wolf et al. (1985). This method requires solvingthe original system and its tangent space (as defined in eq.2) simultaneously with the equations of motion (eq. 1); thatis, the augmented system has n ( n + 1) coupled first-orderequations. There are a number of more recent developmentsthat address some of the numerical difficulties inherent incomputing Lyapunov spectra. The most serious of whichis slow convergence. A strictly positive maximal Lyapunovexponent is often considered as a definition of determinis-tic chaos, but it is difficult to associate them with explicitdynamical quantities (e.g. Cvitanovi´c et al. 2012). More nar-rowly, a positive maximal Lyapunov exponent indicates ex-ponential sensitivity to initial conditions, a feature of anunstable and an irregular system.Despite its wide use, this approach has some intrinsicdisadvantages and interpretive limitations. First, considerthat the general perturbed phase space is foliated by sheaths of chaos around the homoclinic trajectories primary reso-nances with the perturbation (Paper 1, section 4.1 gives asimple example for a toy model). Direct numerical exper-iment has shown that trajectories that are initially in thevicinity of a one of these chaotic sheaths may later on appearregular owing to apparent capture by a regular island; thisstickiness may repeat throughout the computation (Har-soula, Kalapotharakos & Contopoulos 2011). Secondly, ow-ing to the diffusion of the dynamical trajectory through theaction space of the unperturbed phase space, the action val-ues of trajectory in the unperturbed phase space may havedrifted considerably before the time series is long enough tocompute the Lyapunov spectrum. In other words, the long-term converged Lyapunov exponent is telling us somethingabout the evolving dynamical system, but it is hard to relatethis to the finite-time situation typical of galactic dynam-ics. Short-term Lyapunov analysis and the various expan-sion and contraction indices may provide additional usefulinformation, but direct connection to dynamical principlewill remain challenging. Finally, the stochastic layers nearprimary resonances may still result in stochastic behaviourbut remain undetectable by Lyapunov exponents (i.e. weakchaos). See Paper 1 for additional discussion. MNRAS , 1–17 (2015) omputational KAM II These difficulties motivate the exploration of alternativealgorithms for the diagnosis of chaos dynamical systems. Inthe approach to be explored here, we take a constructive ap-proach. Specifically, we begin with the Hamilton-Jacobi (HJ)equation, which is the underlying dynamical framework foraction-angle variables. Paper 1 presents an algorithm for ob-taining a solution to the HJ equation that loosely follows theconstructive proof of the Kolmogorov-Arnold-Moser (KAM)theorem. We interpret the failure to find a solution of theperturbed HJ equation using the KAM construction as in-dicative of a destroyed torus described by some initial value I . This is consistent with the results of explicit integrationpresented in Paper 1. In addition, the numerical KAM (here-after, nKAM) approach will give us expressions for the per-turbed by regular orbits that may then be used for study-ing their implication for galactic structure. For example, ifwe determine the trajectories corresponding to the new al-beit perturbed tori, we can use perturbed tori as basis fornew perturbation theory (see section 5 for additional dis-cussion). Moreover, the nKAM procedure is generally mustfaster than Lyapunov analyses for the same phase-space res-olution. Also, it is embarrassingly parallel (Foster 1995) foran external determined potential field.In some ways, the method to be explored here suffersfrom the same time-scale issue as Lyapunov-exponent com-putation; that is, tori computation for a periodic systemassumes an eternally-applied perturbation, whereas most as-tronomical perturbations have finite duration. It is possibleto use the nKAM method over a finite-time disturbance us-ing extended phase space or employing a Laplace transformin the time domain. This will be explored in a later contri-bution.In this paper, we will apply the nKAM approach to anidealised barred galaxy where the bar is considered to be aperturbation to the axisymmetric galactic disc. We are awarethat this is a tall order for a perturbation theory, but it alsoprovides a good test of its domain of applicability. More-over, barred-galaxy dynamics is a mature subject, largelyexplored with idealised analytic models, Poincar´e surface-of-section plots and Lyapunov spectra. Our hope and goal,here, is the gleaning of some alternative insight from thenKAM approach. Specifically, we will examine the tori bro-ken by adding a quadrupole perturbation to an initially ax-isymmetric exponential disc (see § §
4) as expected (Manos 2008). In the limit of astrong bar (i.e. a large fraction of the disc interior to the baris in the bar), many of the tori inside of the corotation andin bands around ILR and corotation are destroyed, leavinga smaller band of regular orbits of moderate eccentricity leftto support the bar. To be sure, the trajectories for many ofthese broken tori are not dramatically “Brownian”; they ap-pear to flip between different approximate rosettes of vary-ing eccentricity. Presumably, some of these modes for theseweakly chaotic orbits will not be bar supporting all of thetime even if they would be in their unperturbed state. Eventhough the underlying perturbation method may operatingbeyond its domain of applicability, these results suggest that chaos may play a role in limiting the growth and structureof bars, and this is supported by the results of our Lyapunovexponent analysis.
The technical details of the nKAM approach are describedin Paper 1. This section summarises the main ideas andmotivation without the technical and mathematical details.We begin with the standard setting for the equationsof motion in a stellar system: classical Hamiltonian pertur-bation theory (e.g. Binney & Tremaine 2008). The funda-mental coordinate system describing the regular orbits in aquasi-periodic system is the action-angle system. This sys-tem results from a canonical transformation from the 2 n co-ordinates and generalised momenta that separate the gravi-tational potential to those that leave the Hamiltonian a func-tion of n new generalised momenta only. There are n = 3dimensions for galaxies, although we will consider the re-striction to n = 2 for our planar examples here. The equationdefining implicit solution for the generator of the canonicaltransformation that yields the action-angle system is the HJequation (e.g. Goldstein, Poole & Safko 2002). Each regular orbit has a fixed action vector, I , and a corresponding anglevector, w , whose values linearly increase with time. Geo-metrically, the I vector defines a torus in phase space thatbecomes uniformly sampled by the trajectory in time (ex-cept for a small set of degenerate closed-orbit situations).Let us assume that we have a solution of the HJ equa-tion for some background Hamiltonian function that yieldsregular orbits everywhere. Then, let us add a perturbationwith pattern speed Ω p and attempt to find a correction tothe generating function that solves the perturbed HJ equa-tion. If successful, we have a new set of action-angle vari-ables and new tori. We get an algebraic solution for thenew canonical transformation by solving the HJ equationafter a Fourier transform in angle variables. The solutionhas the form of trigonometric series whose arguments havethe form l · w where l is a vector of integers, the Fourierindices. Formally, the coefficients of this series solution has vanishing denominators , indicating the evolution of the spe-cific action values at the point of vanishing is unboundedto linear order. The vanishing denominators have the form (cid:80) j l j Ω j − m Ω p where l j is the Fourier index of the j an-gle variable, m is most often the azimuthal angular har-monic, the Ω j = ∂H/∂I j are the unperturbed frequencies of w which follows from Hamilton’s equations.If each term in the series solution is considered to be in-dependent of the others, the vanishing denominator problemis resolvable. To see this, let us restrict our attention to onlyterm that generates a particular vanishing denominator. Theassociated non-linear problem has a well-defined solution:it is that of a pendulum! The one-dimensional pendulumrepresents the resonant degree of freedom. As individual or-bits exchange angular momentum, the collective changes inactions causes the perturbation to evolve by changing itspattern speed or amplitude, and this drives the perturba-tion the past the resonance for each susceptible trajectory.If we apply the averaging theorem to isolate the evolutionto the specific degree of freedom implied by the vanishingdenominator, and self-consistently include these changes in MNRAS , 1–17 (2015)
M. D. Weinberg the overall solution, we get a prescription for secular evo-lution. Dynamical friction is the classic example of secularevolution of this form (Tremaine & Weinberg 1984). Thistheory is very well developed in celestial mechanics wherethe perturbation strength tends to be very small.Clearly, traditional secular evolution is only piece of thefull dynamical picture. For example, the Kirkwood gaps andsolar system stability cannot be explained using the tools ofsecular evolution. The nature and implications of the full solution remained a mystery until relatively recently in thehistory of astrophysical dynamics. Returning to the previousexample, suppose that we refrain from applying the phaseaveraging to isolate individual terms and retain the wholeseries, vanishing denominators and all. The simple exampleof this is the two interacting resonances considered in Pa-per 1, section 4.1. A series of papers in the early 1960s byKolmogorov, Arnold and Moser (KAM Kolmogorov 1954;Arnold 1963; Moser 1962, 1966) showed that solutions tothe HJ equation nearly always exist as long as we stay farenough away from the phase-space location of the vanishingdenominators. For small perturbations, this is most of phasespace. Note: this does not imply that the vanishing denom-inators are of no consequence (see previous paragraph!) butrather that the dynamics that destroy the tori do not destroythem everywhere in the action-space neighbourhood of theinitial orbit. The general message of KAM theory is: a sig-nificant fraction of tori remain after applying a perturbationas long as the amplitude is sufficiently small.Much of the elegant mathematical machinery for solarsystem dynamics is difficult to apply to problems galacticdynamics for several reasons. First, the perturbations arenot those of point masses but structurally extended withone or more scale lengths. In celestial mechanics, expansionvariables that induce rapidly (i.e. exponentially) converg-ing Fourier series are possible. For example, Laskar (1990)shows that series solutions are practical for planetary per-turbations where the controlling amplitude is the solar massratio and multiplicative terms proportional to powers of ec-centricity and inclination. In galactic dynamics, the pertur-bations are most often globally extended and the series donot converge as quickly. Secondly, the frequency space issimply related to coordinate system in celestial dynamics,that is, the unperturbed frequencies are equal and dependonly on semi-major axis. This not so in galactic dynamicswhere the unperturbed forces are generally non-central andthe frequencies and actions must be computed numerically.Finally, in the decades that followed the publication of theKAM ideas, researchers further elucidated the implication ofthe KAM theory for practical dynamical systems, studyingever larger perturbations and the implications for the struc-ture of real-world phase space. A wide variety of numericalexperiments suggests that the converse of the KAM theoremqualitatively holds: that is, the KAM tori can be destroyedby increasing the strength of the perturbation. This pointis especially relevant to galactic dynamics where perturba-tions in the form of satellite encounters, non-axisymmetricdisturbances (bars and arms) and structural deformations(dark halos and flat discs) are strong perturbations.Our goal here is to use the nKAM approach to identifycircumstances where a large fraction of tori are destroyedand investigate the implication of these broken tori forgalaxy evolution. This study pushes the dynamics beyond that of secular perturbation theory by examining the effectof many terms in the Fourier series solution together. Thatis, we do not isolate and solve each term in the Fourier seriesseparately but include their mutual influence. That is, we al-low may resonances to interact simultaneously. However, tomake a connection with the familiar insights from classi-cal perturbation theory, we will focus on a particular termwhose magnitude is amplified by its proximity to the locus inphase space where the denominator vanishes even though allresonant terms remain in the solution. Following the stan-dard conventions, call this the primary resonance and theremaining multiple resonances secondary resonances. Here,our target application is disc bars considered as a perturba-tion to the axisymmetric equilibrium. There are three strongprimary resonances for a bisymmetric m = 2 harmonic per-turber: 1) the corotation resonance characterised by an orbitmoving at the same azimuthal frequency as the perturbationand inner and outer Lindblad resonances, characterised bythe two radial oscillation for every azimuthal oscillation inthe frame rotating with the bar pattern. The primary reso-nances in the artificial absence of secondary resonances arecharacterised by closed orbits in this rotating frame. Whenisolated to the one-dimensional resonant degree of freedom,this special trajectory that connects the unstable equilibriaand separates the rotation from the libration regions is the homoclinic orbit or separatrix .Owing to the significant strength of bar perturbationsimplied by observations, there are often Fourier coefficientswhose amplitudes are not small relative to the primary res-onance. These secondary resonances may have amplitudesclose to that of the primary resonance, making the distinc-tion somewhat ambiguous. Near the instability for the pri-mary resonance, the secondary terms, even when they arerelatively small, can induce significant changes in the tra-jectory. In effect, the secondary terms noticeably changethe trajectory that would result in absence of these terms.This results in a sheath of stochastic behaviour around theoriginal unstable homoclinic trajectory. For strong pertur-bations, the sheaths of the primary resonance broadenedby secondary resonances may result in large-scale diffusionthrough phase space (i.e. strong chaos, see Zaslavsky 2007).Similar to the overall conclusions of previous work, we do notfind strong chaos to play a major role in barred galaxy dy-namics except, perhaps for nearly radial trajectories and fornear-circular orbits close to corotation. Many authors havestudied this problem in some detail with a variety of novelnumerical techniques (see section 1). For example, Bountis,Manos & Antonopoulos (2012) used a statistical approach toidentify both weakly and strongly chaotic orbits in a barredgalaxy model. Zotos, Caranicolas and collaborators in a se-ries have articles (see e.g. Caranicolas & Zotos 2013) haveused Lyapunov exponents and related indicators to char-acterise chaos in a variety of interesting circumstances. It isnot the goal of this paper to explore the wide range of astro-nomical situations leading to chaos but, rather, to attemptto connect the existence of chaos to standard Hamiltonianperturbation theory and KAM theory.In part, the proof of the KAM theorem rests on therearrangement of terms in the linearised HJ equation toyield a quadratically converging recursive solution for thenew Fourier coefficients. The details of the proof specifysome specific conditions of smoothness and distance in phase MNRAS , 1–17 (2015) omputational KAM II Figure 1.
Circular velocity curve for halo and bulge (blue), thebulge-free halo model (green), the exponential disc (red), the com-bined model including the bulge (cyan), and the combined modelwithout a bulge (magenta). The bulge is constructed to providean approximately flat rotation curve in the inner galaxy. Thebulge-free rotation curve rises inside of approximately two discscale lengths. space from a resonance so that the iteration converges. Theproof also exploits the exponential decrease of coefficient am-plitudes with increasing harmonic order, l . The insight fromKAM is that one may truncate the expansion series at someorder maximal order l max and still recover the correct dy-namics. This is borne out in convergence tests used to verifythe results described below in section 4. Satisfactory con-vergence implies the existence of a torus . In addition, thesolution results from a linearised equation, and it remainspossible that this procedure under- or over-predicts the ex-istence of tori in practice. Finally, it is tempting to interpretthe failure to find a torus as evidence for a ‘broken’ torus,and we will do this here. We have applied the nKAM pro-cedure, the Poincar´e surface-of-section (SOS) method, andLyapunov analysis to the same problem and compared theresults (e.g. section 4 and Paper 1, section 5). Of course,these methods are sensitive to different aspects of the dy-namics, but taken together, the interpretation of ‘broken’tori appears warranted. The numerical details of applyingthis method are described in Paper 1. We adopt the following set of ΛCDM units: the virial radiusis the unit length scale and the mass of the dark matterhalo inside the virial radius is unity . Scaled to the Milky As mentioned previously, our goal is to assess the fraction ofbroken tori, not to show the existence of tori as in the KAM the-orem. To this end, we do not adjust the initial condition in theneighbourhood of the original actions to maintain the Diophan-tine condition that controls the small divisors. In ΛCDM, one uses the uniform spherical collapse model to de-fine radius of collapse, r vir , for region of specified density enhance- Way, the unit mass is approximately 10 M (cid:12) and the unitradius is 300 kpc. We model the bulge and spheroid by twoNFW-like (Navarro, Frenk & White 1997, NFW) profiles. Tocombine the dark-matter halo and bulge models, we “leaveroom” for the inner cusp by adding a core to the dark matterhalo by adopting a modified NFW profile of the form ρ ( r ) = ρ r s ( r + r c ) ( r + r s ) (3)where r c is zero for the original NFW profile. We considera model with and without a bulge. For the model with abulge, the dark-matter halo profile has a concentration of c = 15 ( r s = 0 . M h = 1 inside of r = 1, with r c = 0 .
02. The bulge profile has a concentration of c = 3000( r s = 0 . r c = 1 . ¯6 × − . This coreradius limits the dynamic range of the density inside of anyradius of astronomical interest here which simplifies someof the numerical computations. For the model without abulge, the dark-matter halo profile has a concentration of c = 15 ( r s = 0 . M h = 1 inside of r = 1, with r c = 0 . R ) = M d πa e − R/a (4)The mass of the disc, M d , was adjusted to produce an ap-proximately flat circular velocity curve from zero to ten discscale lengths. This mass is M d = 0 .
04 (approx. 4 × M (cid:12) scaled to the Milky Way). The scale length is chosen to be a = 0 .
01 units (approx. 3 kpc scaled to the Milky Way). Thebulge mass is adjusted to keep the inner profile flat ratherthan rising. This gives a bulge mass 0.02 units (2 × M (cid:12) scaled to the Milky Way) and a B/D ratio of 0.5. The bulge-free model has a rising rotation curve out to several disc scalelengths. The circular velocity profile for each component andthe total for both models is shown in Fig. 1. Motivated byboth observations and the results characterising many pub-lished n-body simulations, we chose the bar length to bethat of disk scale length, a = a . Similarly, we choose pat-tern speed Ω p so that the bar rotates at the same speed asa circular orbit at R = a .Following previous contributions (e.g. Weinberg & Katz2002), we will model the bar as an ellipsoidal mass distri-bution. The gravitational potential for density stratified onellipsoids is well known. Specifically, the gravitational poten-tial internal to the ellipsoid takes the form (Chandrasekhar1969, pg. 52, eq. 93) V = πGa a a (cid:90) ∞ du ∆ (cid:2) Ψ(1) − Ψ( m ( u )) (cid:3) (5) ment. Specifically, the virial mass, M vir , is the mass containedwithin a radius r vir inside of which the mean interior density is∆ times the critical density ρ c : (cid:90) r vir r d r ρ ( r ) = ∆3 ρ c r where ρ ( < r vir ) = ∆ ρ c is the halo’s average density within thatradius. Cosmological n-body simulations are used to calibrate ∆.The virial radius may be related to point within which the ma-terial obeys the virial relation and external to which the mass isstill collapsing onto the object.MNRAS , 1–17 (2015) M. D. Weinberg where m ( u ) = (cid:88) i =1 x i a i + u (6)and Ψ( m ) = (cid:90) m dm ρ ( m ) . (7)The gravitational potential external to the ellipsoid has theform (Chandrasekhar 1969, pg. 51, eq. 89) V = πGa a a (cid:90) ∞ λ du ∆ (cid:2) Ψ(1) − Ψ( m ( u )) (cid:3) . (8)We considered three different types of densities, powerlaw, Ferrers (1887), and exponential: ρ ( m ) = ρ (cid:40)(cid:0) − log m (cid:1) if γ = − (cid:0) − m γ (cid:1) otherwise ρ ( m ) = ρ (cid:0) − m (cid:1) γ ρ ( m ) = ρ (cid:16) e − /γ − e m /γ (cid:17) with γ chosen to be physically sensible (i.e. γ > = 0 for thefirst two cases and γ of order the disc scale length in thefinal case). The resulting quadrupole was similar with allthree models, and for ease, we choose to represent the barwith a homogeneous ellipsoid with axes a = a = 0 . , a =0 . , a = 0 . U ( r ) = b r (cid:16) rb (cid:17) (9)Fig. 2 compares the fit from equation (9) to the quadrupolecomponent of equations (5) and (8). The power-law fit sys-tematically exceeds the true value for small r (cid:28) .
01 butotherwise captures the run of U in the vicinity of the barend at 0.01 quite well. In this section, we use the nKAM method developed in Pa-per 1 and the bar perturbation described in section 3 todetermine the location and examine the morphology of bro-ken tori. The approach is that same as that described forthe Kepler example in Paper 1, section 5. The unperturbedmodel is the entire axisymmetric component: bulge, dark-matter halo and the exponential disc. The quadrupole fit tothe triaxial ellipsoid representing the bar, described in theprevious section, is the perturbation. Thus, the perturbationamplitude only affects the quadrupole strength. We assumea constant pattern speed Ω p fixed to that of the circularorbit frequency at one disc scale length. As mentioned insection 2, it is possible to use the nKAM method to exploreperturbations of finite duration; these will be tackled in alater contribution.Since the actions are the fundamental description of theunperturbed phase space, we describe results of the nKAM -6 -5 -4 -3 -2 -1 P o t e n t i a l exactfit -3 -2 -1 radius−0.03−0.02−0.010.000.010.020.03 D i ff e r e n c e Figure 2.
Comparison of the exact quadrupole potential U ( R )from the homogeneous ellipsoid and the power-law-type fit (eq. 9.The upper panel shows compares the fit in logarithmic units andlower panel shows the difference between the fit and the exactquadrupole. procedure in the radial action–azimuthal action plane (thatis, the I – I plane) for a variety of bar quadrupole ampli-tudes. The fiducial mass in the triaxial ellipsoid that repre-sents the bar is set equal disc mass enclosed within the barradius. The resulting quadrupole amplitude is scaled to thedesired bar strength. For visual comparison, Fig. 3 shows themidplane density implied by the quadrupole perturbationfor four bar strengths spanning the range of interest. We usea logarithmic scaling for comparison with astronomical im-ages and to better illustrate the bar amplitude near the endsof the bar. The full-amplitude ( (cid:15) = 1) profile bears a close re-semblance observed luminosity density for strong bars (e.g.Gadotti & de Souza 2006). The negative-density “dimples”perpendicular to the bar that occurs for strongest amplitudecase is an artefacts of only including the quadrupole term;higher-order multipole terms would fill in these regions per-pendicular to the bar’s major axis. The potential and forceremain physical even when dimple artefacts appear in thedensity profile. Fig. 3 further suggests that perturbationswith (cid:15) (cid:38) . (cid:15) = 0 . , . , . , .
0. Each filled point denotes a broken torus; the colour key labels eachby the term in the perturbation series (see Paper 1, eq. 16)with maximum amplitude. The loci of broken tori follow thelow-order resonances, as expected. For example, the locus ofbroken tori at the lower-left in Fig. 6b follows the corotationresonance (0 , , MNRAS , 1–17 (2015) omputational KAM II (a) (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 1 . Figure 3.
Midplane density of the exponential disc with scale length a = 0 .
01 (system units) including the bar perturbation for severalamplitudes, (cid:15) . The bar scale length equals the disc scale length, a = 0 . ten the inner Lindblad resonance (ILR, l = − , l = 2),and vice versa for the primary ILR. For (cid:15) = 0 . (cid:36) l of the W series fall between 0.3% and 3%; these amplitudes seemto increase roughly linearly with increasing (cid:15) . As the am-plitude continues to increase, the width of the region nearthe locus of corotation, OLR, and ILR increase, starting toclose the gap of unbroken tori between ILR and corotationavailable to support the bar figure. The main differences be-tween the two series, Figs. 6a–6d and Figs. 7a–7d, are thelocations and existence of various resonances. A flat rotationcurve results in more resonances near and around vicinity ofthe bar, including the ILR at small radii.We will see in section 4 below that the nKAM compu-tation appears to underestimate the measure of stochasticorbits for mild eccentricity, typical of a disc, suggesting thata large fraction of orbits in the bar region are stochastic atsome level. Indeed, computation of the SOS plots for in thisregion suggests most of the orbits in this region with near circular orbits, I (cid:28) I , show some stochastic spreading.In many but not all cases, these irregular orbits appear tohave inner and outer turning points although the envelopeof the trajectory vary with time. This calculation is not selfconsistent so the we have no way of assessing whether theremaining unbroken tori are capable of reproducing the im-posed gravitational perturbation. Nonetheless, the nKAMresults suggest that a large fraction of the phase space inthe bar vicinity is chaotic owing to bar itself. This raisesthe possibility that the bar amplitude and structure may beself-limited by stochasticity.We now compare the results of a Lyapunov exponentanalysis with nKAM predictions. The presence of a bulgeadds an additional interpretive complication so we begin bydescribing the results for the bulge-free galaxy model with arising rotation curve. The regions of positive Lyapunov expo-nent for the amplitudes in Figs. 6a–6d are depicted in Figs.8a–8d. The locus of positive exponents and broken tori coin-cide approximately for the corotation and OLR resonances MNRAS , 1–17 (2015)
M. D. Weinberg
Figure 4.
Mapping of action space to orbit turning points for the unperturbed model shown on large (left) and small (right) scale. Lociof inner turning points (blue) and outer turning points (red) are labelled in system units. For reference, the bar length and disc scalelength is 0.01 system units.
Figure 5.
Twenty consecutive orbits from the two-resonancemodel from Paper 1 whose surface of section is illustrated section4.1, Fig. 2 for an initial condition close to the homoclinic trajec-tory of the primary resonance. As the trajectory approaches theunstable equilibrium of the primary resonance at ( I , w ) = (0 , at all amplitudes. As the bar amplitude increases, corota-tion appears in positive Lyapunov values at high ( I (cid:38) I )and low ( I (cid:28) I ) eccentricity. Although the ILR does notappear as a primary resonance, it is a secondary resonancenear the inner edge of the region of broken tori correspond-ing to the corotation resonance in all of Figs. 6a–6d. At loweramplitude, the ILR shows up as a distinct locus (Figs. 8aand 8b). At higher amplitude, regions of positive Lyapunovexponent owing to corotation and ILR blend together in asingle locus (Figs. 8c and 8d).For the flat rotation curve model, the regions of positiveLyapunov exponent for the amplitudes in Figs. 7a–7d are de- picted in Figs. 9a–9d. For all amplitudes, we find a qualita-tive demarcation at I ≈ . ∝ r , stronglyaffects the bulge dynamics. The perturbation causes the lo-cation of these high-frequency striations to change locationin action space with bar strength. In essence, the striationsmanifest in Figs. 9a–9d are a form of orbit “shocking.”The locus of positive exponents and broken tori coincideapproximately for the ILR at all amplitudes. SOS plots con-firm that these orbits exhibit significant irregularity. As thebar amplitude increases, corotation appears in positive Lya-punov values at high ( I (cid:38) I ) and low ( I (cid:28) I ) eccentric-ity. All in all, the overall conclusions for both galaxy modelsare qualitatively the same: increasing amplitude gives riseto an increasing fraction of stochastic orbits in the bar re-gion. The galactic model with a flat rotation curve has alarger zone of unbroken tori that are bar supporting thanthe model with a rising rotation curve. For both models,the nKAM analyses predict a region at moderate eccentric-ity ( I (cid:46) I ) just below the chaotic corotation regions wherethe orbits remain approximately regular in morphology andmay support the bar. This zone of regularity decreases withincreasing bar amplitude.A major advantage of the nKAM approach is that itnaturally provides a dynamical characterisation for the bro-ken torus, i.e. hints to the primary and secondary reso-nances involved. On the other hand, a broken torus pre-dicted by the KAM method does not predict the degreeof exponential sensitivity. However, empirically, a ratio ofsecondary to primary amplitude of the generating functioncoefficients (cid:36) l (Paper 1, eq. 16) of order 0.1 or larger nearly MNRAS , 1–17 (2015) omputational KAM II (a) (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 1 . Figure 6.
Broken (good) tori are indicated by filled circles (dots) as a function of the unperturbed values of radial action I andazimuthal action I for the bulge-free galaxy model. The broken tori are colour coded by the principle resonance in the form ( l , l )(legend). For example ( − ,
2) is the inner-Lindblad resonance and (0 ,
2) is the corotation resonance. Initially circular orbits like alongthe I = 0 vertical line and initially radial orbits like along the I = 0 horizontal line. The right-hand vertical axis describes the radiusof the circular orbit ( I = 0) corresponding to each value of I . always has a numerically measurable positive Lyapunov ex-ponent. In essence, the secondary resonances break the for-ward and backward time symmetry of the stable and un-stable branches of the resonant orbit. Chaos is induced byrapidly oscillating intersections of the perturbed stable andunstable manifolds in the sense of the Birkhoff-Smale theo-rem (Smale 1965) that gives rise to exponential sensitivity.The resulting chaos is confined to a sheath surrounding theoriginal unperturbed homoclinic trajectory (as in Holmes &Marsden 1982). The dynamics of this behaviour was con-sidered in detail by Wiggens (1990, Chapter 4). It can beeasily seen numerically by plotting the solution of a directnumerical integration for the many resonance system in theaction-angle space of the one-dimensional primary resonancenear the homoclinic trajectory. For example, Fig. 5 shows aThe resulting motion for broken tori of this type re-mains generally confined in phase space. However, as the perturbation grows, the chaotic sheath can fill a large frac-tion of the original libration zone, which may be large. Atypical orbit in this sheath appears to quickly vary its eccen-tricity; visually it presents as a rosette with turning pointsvarying in radius with time.This insight helps us interpret the differences betweenthe approaches. Recall that the Lyapunov value indicatesthe exponential rate of divergence of two initial conditions.However, this does not describe the behaviour of the tra-jectory. Qualitative classification of orbital regularity by eyecombined with the computation of SOS plots suggests signif-icant irregularity for exponents with λ = O (10 − ). In manycases, the nKAM locus provides a more sensitive indicatorof separatrix disruption than the Lyapunov exponent. Forexample, Fig. 10 depicts the trajectories and SOS plots fortwo orbits with the same value of I and nearly the samevalue of Lyapunov exponent. One orbit is in the chaotic zone MNRAS , 1–17 (2015) M. D. Weinberg (a) (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 1 . Figure 7.
As in Fig. 6 but for the galaxy model with a bulge and flat rotation curve. parented by corotation as predicted by the nKAM method,and one orbit is between the ILR and corotation chaoticzones [cf. Fig. 7b]. The orbit that is predicted to be regularby nKAM shows qualitatively regular features, and this iscorroborated by the SOS plot. The orbit that is predicted tobe chaotic by nKAM has two modes; the perturbation hasbroken the original torus resulting into two loop-like tra-jectories that joined at an unstable point. This behaviouris typical of broken tori in the corotation loci of Fig. 7 for I ∼ > I . Broken tori with I (cid:28) I exhibit more complex be-haviour, presumably from the superposition of effects frommultiple secondary terms. This “mode switching” explainsthe qualitatively ‘chaotic’ appearance of trajectory but thesmall measured Lyapunov exponent values. Only near thehyperbolic point of the islands in the stochastic sheath doesthe exponential sensitivity to phase space play a strong role,and therefore, the time-averaged Lyapunov exponent can re-main anomalously small in light of the qualitative stochasticbehaviour of the trajectory over many periods.In summary, each of these two methods, nKAM andLyapunov exponent computation, reveal different aspects of irregular systems. The nKAM procedure identifies locationsin phase space where influence of other nearby or strong res-onances are capable of inducing morphological changes innear resonant orbits. Lyapunov-exponent analysis attemptsto measure the rate of exponential sensitivity directly. Thistechnique is good at discovering strong chaos but may misssubtle cases where the mild break up of a resonance into sev-eral islands that permits an orbit to spontaneously changemorphologies after many dynamical times. In this later case,the trajectory is only exponentially sensitive to its initialconditions near an unstable point, and the time to trans-verse this critical region is a small fraction of its orbital time.While in this critical region, the orbit may change from onenearly regular morphology to another. In other words, thecharacteristic time scale for exponentiation, measured usingthe standard Lyapunov exponent algorithm, may be muchlonger than an a Hubble time but the broken torus can havea qualitatively different character owing to its tendency toswitch between morphologies.The nKAM method does appear to underpredict thestochastic regions. Indeed, many of the orbits with I < MNRAS , 1–17 (2015) omputational KAM II (a) (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 1 . Figure 8.
Values of the largest (positive) Lyapunov exponents for models with various perturbation amplitudes for the bulge-free galaxymodel. Positive Lyapunov exponents are determined by regression analysis to determine convergence with increasing time as described inPaper 1, section 3.2. The local estimate is computed by smoothing the ensemble onto a regular grid. Contour values are logarithmicallyspaced; values of ( I , I ) that are converging to zero are shown in white. .
001 and 0 . < I < .
015 that show larger values ofLyapunov exponent in Fig. 9 exhibit the qualitative signa-ture of chaos, although the nKAM method is able to con-struct a torus. Choosing a specific example from the (cid:15) = 0 . I , I ) = (0 . , . conservative : itwill find regular orbits if they are in a small neighbourhoodof the specified initial condition even if irregular orbits can be found there as well. Conversely, the Lyapunov exponentprocedure seems to find both regular and irregular orbitsin the same neighbourhood, if the both exist, and possiblyswitches between them (more on this below).A second limitation of the nKAM procedure using thefundamentally linearised solution of the HJ equation is its in-trinsically perturbative nature. So, it is possible that nKAMmethod has been pushed beyond its domain of accuracy forthis problem resulting in the underprediction of stochastic-ity in this case. For larger value of I , the nKAM appearsto be a sensitive indicator of weak or confined (e.g. Dvo-rak et al. 1998; Winter, Mour˜ao & Giuliatti Winter 2010)chaos, perhaps as expected, than the Lyapunov exponentvalue. However, the different characterisations persist evenfor (cid:15) (cid:46) .
2, and therefore, so strong non-linearity is not thefull explanation.For another example, Fig. 12 shows the same plots asFig. 10 for a chaotic orbit in the ILR region. Upon increasing
MNRAS , 1–17 (2015) M. D. Weinberg (a) (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 1 . Figure 9.
As in Fig. 8 but for the galaxy model with a bulge and flat rotation curve. I at fixed I beyond the locus of nKAM-predicted brokentori, the SOS continues to exhibit significant width, suggest-ing that the nKAM method underestimates the measure ofbroken tori here as well. Decreasing the amplitude of the per-turbation reveals the bifurcations and islands expected of astochastic layer; the trajectory appears regular over shorttimes but can switch between different regular morpholo-gies. This reinforces our previous interpretation for the dis-crepancy between the apparently small values of Lyapunovexponent when the orbit appears irregular overall: the stick-iness in the stochastic sheath results in regular behaviourpunctuated by jumps to different parent resonances near is-lands. One may easily verify there is clear structure in theSOS over short finite-time periods that blends together overlong periods. In this paper, we using the numerical KAM (nKAM) methoddescribed in Weinberg (2015a) to explore the location andfrequency of tori broken by a Galactic-bar perturbation. Inessence, this technique predicts the existence of regular or-bits (which have three actions) and irregular orbits (withfewer than three actions) by an iterative solution of the lin-earised Hamilton-Jacobi (HJ) equation. The HJ equationdefines a canonical transformation to the action-angle vari-able system, and we assume that a consistent solution forthe canonical formation implies existence and vice versa. Insection 4, we show that the nKAM computation yields usefulpredictions predictions for regions of irregularity, however,the method underpredicts the width of the chaotic zonesfor strong bars. This trend may be understood as follows.The nKAM method appears to find a torus if one exists inthe neighbourhood of the initial conditions, even if brokentori exist in this same neighbourhood. Conversely, the Lya-punov method tends to find exponential divergence in the
MNRAS , 1–17 (2015) omputational KAM II (a) regular zone, trajectory (b) regular zone, sos(c) chaotic zone, trajectory (d) chaotic zone, sos Figure 10.
Orbits shown in the frame rotating with the bar for (cid:15) = 0 .
3. The regular orbit constructed by the nKAM method has( I , I ) = (0 . , . I , I ) = (0 . , . λ ≈ × − . The left-hand panels plot a short segment of the trajectory used to compute the two-sided SOS plots in the right-handpanels. same circumstances, if it exists, although orbits may switchbehaviour suddenly, e.g. if the initially diffusing trajectoryencounters a regular island. These differences are not ex-pected, but further elucidate that chaos is a broad catch-allterm for orbits with a wide variety of irregular morphology.Thus, the two approaches reveal different informationabout the same problem. The nKAM provides dynamicaldetails about the principle and secondary resonances forbroken tori although the method is less accurate for strongperturbations owing to its perturbative nature. The Lya-punov analysis works for any perturbation, but the dynam-ical origin of a positive Lyapunov exponent is not readilyavailable. Moreover, our investigation suggests that Lya-punov exponent analysis does not provide useful informationabout stochastic orbits that appear isolated to stochasticlayers. The morphology of the orbits within these layers of-ten exhibit large variance in orbital frequencies while yield-ing the anomalously small Lyapunov exponents typical of weak chaos. Since the morphology of orbits as measured bytheir time-averaged density is key to understanding the evo-lution of patterns such as arms and bars, an random switch-ing between families matters. This suggests that Lyapunovanalysis may have less value for some important galactic dy-namics problems. In addition, the nKAM is capable of pro-viding a new action-angle coordinate system for those torithat survive a perturbation, although we do not exploit thisfeature in the present paper. Conversely, the nKAM methoddoes not provide a direct assessment of the magnitude of thechaos without additional analysis; e.g. strong chaos may beinferred from the relative magnitude of coefficients in (cid:36) l using the Chirikov (1979) criterion.Our conclusion from this work is that galactic evolutionunder the influence of strong patterns (bars, spiral arms,etc.) has an important chaotic component. As described insection 1, this conclusion has been reached by previous re-searchers as well. Although standard perturbation theory MNRAS , 1–17 (2015) M. D. Weinberg
Figure 11.
Orbit (left) and two-sided Poincar´e SOS plot (right) for an initially nearly circular orbit with ( I , I ) = (0 . , . (cid:15) = 0 . Figure 12.
Orbit (left) and two-sided Poincar´e SOS plot (right) for an initially mildly eccentric orbit with ( I , I ) = (0 . , . does remain isolated. provides some physical insight, it is necessarily an incom-plete picture of the underlying dynamics. We have somehope that the extended perturbation theory described abovemay provide a bridge between analytic methods and un-encumbered n-body simulation. Numerically, much of ourmotivation for criteria for orbit integration stems from thestudy of regular orbits. In addition to characterising the new method, our maingoal is a study of the degree of stochasticity induced bya stellar bar in an initially axisymmetric galaxy disc. Theimportance of weak chaos in structuring stable orbits un-der non-axisymmetric perturbations and the underlying dy-namics is nicely illustrated by Romero-G´omez et al. (2006,2007); Athanassoula, Romero-G´omez & Masdemont (2009)who describe the importance of weak chaos in creating ringsin barred galaxies. Their dynamical description of orbit tran- sitions at unstable points is essentially that described insection 4 (Fig. 10d and associated discussion.). There are,of course, other approaches. For example, Bountis, Manos& Antonopoulos (2012) describe a probabilistic approachthat simultaneously discriminates between weak and strongchaos. An advantage of the nKAM approach is that it pro-vides a numerical tool for a wholesale investigation of phasespace as a in a computationally efficient way. At the sametime, it produces a multiple-resonance model Hamiltonianto further investigate broken tori.Based on an action-space expansion, our nKAM inves-tigation provides a provides a statistical view of irregularitycaused by a perturbation or pattern. We assumed an ana-lytic bar model and examined the influence of its quadrupolemoment on the disc. The bar here is modelled by a homoge-neous ellipsoid with the axis ratio and density profile typicalof observed bars (see section 3). The bar does not evolve withtime or is it self-consistent in any way; this is an idealisationthat makes this present investigation tractable, although ex-
MNRAS , 1–17 (2015) omputational KAM II tensions with more generality are possible. The key findingsare as follows:(i) The principle resonances causing broken tori in theinner galaxy are the corotation (0 ,
2) and the OLR (1 , − ,
2) and OLR (1 , , As described by Chirikov (1979), in Paper 1, section 4.1,and in the previous sections of this paper, non-deterministicbehaviour owing to mutual interaction of resonances givesrise to a sheath of chaotic behaviour around the primaryresonance. In perturbation theory based on the averagingprinciple, we isolate the degree of freedom correspondingto a resonance by examining its dynamics close to its (un-stable) homoclinic trajectory. The phase frequency corre-sponding to this degree of freedom will be small and wecall the corresponding action the slow action. The slow ac-tion corresponding to this homoclinic trajectory of the pri-mary resonance will be, in general, a linear combination ofthe actions defined by the separable coordinates of the un-perturbed background potential. Therefore, motion in thissheath can result in changes of the original actions of orderthe width of the libration zone.For the specific example explored in section 4, chaot-ically induced switching of orbit families occurs in regionsof very small phase-space volume near unstable points. Thepervasiveness of this fine-scale structure in phase space in-duced by chaos motivates a technical campaign to determinesufficient conditions on n-body simulation necessary to ac-curately recover the important consequences of the mixedregular and irregular dynamics of galaxies. Typical time-step selection algorithms for n-body simulations choose afraction of the rate of change of position or velocity of atrajectory in its gravitational field. The scale of the gravita-tional field may be determined by examining the derivativeof forces or from prior knowledge about the desired reso-lution length scale. However, the dynamics of interacting resonances is sub-dimensional. Therefore, the details of thedynamics around the primary resonance may require the res-olution of scales much smaller than those characteristic ofthe trajectory as a whole. For example, if a single step in theODE solver moves the orbit through the chaotic zone alto-gether, the chaotic dynamics may disappear from the simu-lation. Indeed, accurate computation of the Lyapunov expo-nents in Paper 1, section 3.2 suggest that the required timesteps are much smaller than those chosen for n-body sim-ulations; larger time steps resulted in negative rather thanpositive exponents in many cases. Similar issues were seen inthe computation of SOS plots for section 4. The topology ofthe SOS plots for broken tori are often not converged untiltime steps are several orders of magnitude smaller than thecharacteristic orbital times for the RK4 time stepper.
Many of our analytic tools for studying modes and sec-ular processes excited by bars and strong spiral patternsare based on the existence of regular orbits. Our success instudying the statistical implication of an imposed patternas described in section 5.2 leads to a secondary longer-termgoal: can we extend perturbation theory methods to under-stand and predict the long-term future of secularly evolvinggalaxies that includes the effects of chaos?Indeed, an immediate consequence of the unbroken torifound using nKAM is a new basis for secular perturbationtheory. The framework described in Paper 1 and in section2 provides the new regular trajectories in an action-anglerepresentation when they exist. This basis may be used toestimate bar kinematic signatures, perform stability analy-ses, and so on; any task previously performed using Hamil-tonian perturbation theory may be performed using the newaction-angle representation. This naturally allows the stan-dard tools of secular perturbation theory (i.e. Tremaine &Weinberg 1984) to be used in the mixed context of regu-lar and irregular orbits. Of course, this only makes sense ifmany of the tori remain unbroken. Based on the results ofsection 4, this roughly true for all but the strongest bars.The orbit morphology in the perturbed regions may besignificantly different than the unperturbed one. For exam-ple, Fig. 10a describes an regular orbit reconstructed by thenKAM method. The SOS plot [Fig. 10d] reveals that thesimple rosette orbit has been destroyed and replaced withan orbit whose turning points are trapped by the rotat-ing bar potential. These changed morphologies of nKAM-determined regular orbits are capable of supporting newclasses of responses not possible with the original unper-turbed regular system. For Fig. 10c, the new orbit may stillbe bar supporting, but less so than a classical x1 bar orbitfrom Contopoulos & Papayannopoulos (1980). In addition,each new class of regular orbits is likely to present a differ-ent kinematic signature than the original unperturbed orbit,and an ensemble chosen to represent the bar may be usedto predict kinematic signatures.A companion paper applies the nKAM method to athree-dimensional but axisymmetric system. A natural ex-tension of this work to three dimensions would explore therole stochasticity in coupling between the planar and ver-tical degrees of freedom. For some examples of the ques-tions to be addressed, does the bar perturbation drive disc
MNRAS , 1–17 (2015) M. D. Weinberg thickening through stochastic means? Do broken tori playa role in producing a pseudo-bulge? As described in Pa-per 3, the fully asymmetric three-dimensional implementa-tion of the nKAM method is computational challenging butsurmountable. Each torus is represented by several three-dimensional Fourier series describing the canonical generat-ing function and tables of the unperturbed orbits as a func-tion of angles. The numerical partial derivatives requiredby the nKAM method requires that three-dimensional ac-tion cube be finely sampled by these tori. This leads to alarge RAM requirement. The current code parallelises toruscomputation and shares the results with all processes. A re-designed code would use an action-space domain decompo-sition to eliminate need for each compute node have a fullcopy of the Fourier series and angles grids for each torus.This will allow a computationally feasible, general, three-dimensional nKAM code.The example here assumes the simplest quasi-periodictime-dependent perturbation: one with a single constantpattern speed. As described in Paper 1 and sections above,the nKAM method may be extended to include generaltime-dependent perturbations by breaking down the time-dependent perturbation using a Fourier (or Laplace) seriesapproximation. A related question is the effect of the slowevolution of the bar and background gravitational field ow-ing to secular torques. For example, does the slow evolutionincrease or decrease the size and effects of the stochasticlayers? Naively, the effective frequency of the slow evolutionis much smaller than the orbital frequencies in the innergalaxy, and this makes the coupling between the naturalfrequencies and the evolution frequencies weak. Conversely,the low spatial frequency of the secular evolution may ef-fectively couple to the primary perturbation for a strongbar, providing new channels for separatrix breaking. More-over, the existence of chaos in the vicinity of orbits sup-porting the bar is likely to affect capture into libration. Forexample, since bars are formed by capture into libration, wemight expect diffusion in and out of the libration zone as thebar pattern speed changes. Results from section 4 suggestthat stochastic regions caused by the forming and growingbar increase with bar strength. Understanding influence ofslow evolution on stochasticity from first principles and thepossibility that bars are chaotically self-limited is the nextchallenge.
ACKNOWLEDGMENTS
This work was supported in part by NSF awards AST-0907951 and AST-1009652. MDW gratefully acknowledgessupport from the Institute for Advanced Study, wherework on this project began. Many thanks Douglas Heggiefor many valuable comments on an early version of thismanuscript.
REFERENCES
Arnold V. I., 1963, Russian Mathematical Surveys, 18, 13Athanassoula E., Romero-G´omez M., Bosma A., Masdemont J.,2009, Monthly Notices of the Royal Astronomical Society, 400,1706 Athanassoula E., Romero-G´omez M., Masdemont J., 2009,Monthly Notices of the Royal Astronomical Society, 394, 67Benettin G., Galgani L., Giorgilli A., Strelcyn J. M., 1980, Mec-canica, 15, 9Binney J., Tremaine S., 2008, Galactic Dynamics: (Second Edi-tion) (Princeton Series in Astrophysics). Princeton UniversityPressBountis T., Manos T., Antonopoulos C., 2012, Celestial Mechan-ics and Dynamical Astronomy, 113, 63Brunetti M., Chiappini C., Pfenniger D., 2011, A&A, 534, A75Cachucho F., Cincotta P. M., Ferraz-Mello S., 2010, Celestial Me-chanics and Dynamical Astronomy, 108, 35Caranicolas N. D., Zotos E. E., 2013, Publ. Astron. Soc. Australia,30, 49Chandrasekhar S., 1969, Ellipsoidal Figures of Equilibrium. YaleUniversity Press, New HavenChirikov B. V., 1979, Physics Reports, 52, 265Contopoulos G., 2002, Order and chaos in dynamical astronomy.Order and chaos in dynamical astronomy /George Contopou-los. Berlin : Springer, c2002. (Astronomy and astrophys li-brary) QB 351 .K58 2002Contopoulos G., Harsoula M., 2010, Celestial Mechanics and Dy-namical Astronomy, 107, 77Contopoulos G., Papayannopoulos T., 1980, A&A, 92, 33Cvitanovi´c P., Artuso R., Mainieri R., Tanner G., Vattay G.,2012, in Chaos: Classical and Quantum, verstion 14 edn., NielsBohr Institute, Copenhagen, pp. 114–124Dvorak R., Contopoulos G., Efthymiopoulos C., Voglis N., 1998,Planet. Space Sci., 46, 1567Ferrers N. M., 1887, Quart. J. Pure Appl. Math, 14, 1Foster I., 1995, Designing and Building Parallel Programs: Con-cepts and Tools for Parallel Software Engineering. Addison-WesleyGadotti D. A., de Souza R. E., 2006, ApJS, 163, 270Giordano C. M., Cincotta P. M., 2004, A&A, 423, 745Goldstein H., Poole C. P., Safko J., 2002, Classical Mechanics,3rd edn. Addison-WesleyHarsoula M., Kalapotharakos C., Contopoulos G., 2011, MNRAS,411, 1111Holmes P., Marsden J. E., 1982, Comm. Math. Phys. 82 523(1982), 82, 523Kaufmann D. E., Contopoulos G., 1996, A&A, 309, 381Kolmogorov A., 1954, The Proceedings of the USSR Academy ofSciences, 98, 527530Laskar J., 1990, Icarus, 88, 266Manos A., 2008, PhD thesis, Universit/’e de Provence, Franceand University of Patras, GreeceManos T., Athanassoula E., 2011, Monthly Notices of the RoyalAstronomical Society, 415, 629Manos T., Machado R. E. G., 2014, MNRAS, 438, 2201Moser J., 1962, Nachrichten der Akademie der Wissenschaften inGttingen, Mathematisch-Physikalische Klasse, 1Moser J., 1966, Annali della Scuola Normale Superiore di Pisa,Ser. III, 20, 499Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Patsis P. A., Athanassoula E., Quillen A. C., 1997, ApJ, 483, 731Romero-G´omez M., Athanassoula E., Masdemont J. J., Garc´ıa-G´omez C., 2007, A&A, 472, 63Romero-G´omez M., Masdemont J. J., Athanassoula E., Garc´ıa-G´omez C., 2006, A&A, 453, 39Shimada I., Nagashima T., 1979, Prog. Theor. Phys., 61, 1605Smale S., 1965, in Differential and Combinatorial Topology (ASymposium in Honor of Marston Morse), Princeton Univer-sity Press, Princeton, N.J., p. 6380Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729Tsoutsis P., Kalapotharakos C., Efthymiopoulos C., ContopoulosG., 2009, A&A, 495, 743Weinberg M. D., 2015a, MNRAS, paper 1, submittedMNRAS , 1–17 (2015) omputational KAM II Weinberg M. D., 2015b, MNRAS, paper 3, submittedWeinberg M. D., Katz N., 2002, ApJ, 580, 627Wiggens S., 1990, Introduction to Applied Nonlinear DynamicalSystems and Chaos, Texts in Applied Mathematics. SpringerWinter O. C., Mour˜ao D. C., Giuliatti Winter S. M., 2010, A&A,523, A67Wolf A., Swift J. B., Swinney H. L., Vastano J. A., 1985, Physica,16D, 285Zaslavsky G. M., 2007, The Physics of Chaos in HamiltonianSystems. Imperial College PressMNRAS000