Manifold spirals and secular evolution in galaxies

Manifold spirals, disc-halo interactions and the secular evolution in N-body models of barred galaxies

C. Efthymiopoulos P.E. Kyziropoulos, R.I. Páez K.Zouloumi and G.A. Gravvanis
Research Center for Astronomy, Academy of Athens, Soranou Efessiou 4, 11527 Athens, Greece
Department of Electrical and Computer Engineering, School of Engineering,
Democritus University of Thrace, University Campus, Kimmeria, 67100 Xanthi, Greece
Department of Mathematics, University of Padova, Via Trieste 63, 35121 Padova, Italy
Department of Physics, University of Athens, Panepistimiopolis, 11521 Athens, Greece
E-mail: cefthim@academyofathens.grE-mail: pkyzirop@ee.duth.grE-mail: paez@math.unipd.itE-mail: konstantina-z7@hotmail.comE-mail:
Accepted XXX. Received YYY; in original form ZZZ

The manifold theory of barred-spiral structure provides a dynamical mechanism explaining how spiral arms beyond the ends of galactic bars can be supported by chaotic flows extending beyond the bar’s co-rotation zone. We discuss its applicability to N-body simulations of secularly evolving barred galaxies. In these simulations, we observe consecutive ‘incidents’ of spiral activity, leading to a time-varying disc morphology. Besides disc self-excitations, we provide evidence of a newly noted excitation mechanism related to the ‘off-centering’ effect: particles ejected in elongated orbits at major incidents cause the disc center-of-mass to recoil and be set in a wobble-type orbit with respect to the halo center of mass. The time-dependent perturbation on the disc by the above mechanism correlates with the excitation of new incidents of non-axisymmetric activity beyond the bar. At every new excitation, the manifolds act as dynamical avenues attracting particles which are directed far from corotation along chaotic orbits. The fact that the manifolds evolve morphologically in time, due to varying non-axisymmetric perturbations, allows to reconcile manifolds with the presence of multiple patterns and frequencies in the disc. We find a time-oscillating pattern speed profile at distances between the bar’s corotation, at resonance with the succession of minima and maxima of the non-axisymmetric activity beyond the bar. Finally, we discuss disc thermalization, i.e., the evolution of the disc velocity dispersion profile and its connection with disc responsiveness to manifold spirals.

galaxies: structure – kinematics and dynamics – spiral
pubyear: 2019pagerange: Manifold spirals, disc-halo interactions and the secular evolution in N-body models of barred galaxiesA

1 Introduction

The manifold theory of spiral structure (Romero-Gomez et al. (2006); Voglis et al. (2006a); see also Danby (1965)) provides a dynamical mechanism to interpret how structures such as the bi-symmetrical spiral arms, observed beyond the ends of a galactic bar, can arise out of chaotic motions of the underlying distribution of matter. Two versions of the manifold theory have been developed so far in the literature.

i) In the ‘flux-tube’ version (Romero-Gomez et al. (2006); Romero-Gomez et al. (2007); Athanassoula et al. (2009a); Athanassoula et al. (2009b); Athanassoula et al. (2010); Athanassoula (2012)), one considers continuous-in-time orbits describing the flow of matter away from the bar’s unstable Lagrangian points L and L, as viewed in a frame of reference rotating with angular speed equal to the bar’s pattern speed . The unstable ‘flux-tube’ manifolds are invariant sets formed by all orbits tending asymptotically towards or in the backward sense of time. This means that, in the forward sense of time, these orbits form an outflow directed away from or . An elementary linearization of the equations of motion around or shows that these outflows, or ‘flux-tubes’, have the form of trailing spiral arms. Besides and , similar flux-tube manifolds can be constructed for the whole family of epicyclic periodic orbits (called ‘Lyapunov orbits’) around or .

ii) In the version of the manifold theory called ‘apocentric manifolds’ (Voglis et al. (2006a); Tsoutsis et al. (2008); Tsoutsis et al. (2009); Efthymiopoulos (2010); Harsoula et al. (2016)) one computes first the flux-tube manifolds, and then isolates only those points which correspond to apsidal positions, i.e., local apocentric or pericentric points of each orbit in the flux-tube. One can show that the pattern formed by the union of the apocentric points yields again trailing spiral arms. Near or the shapes of the flux-tube and apocentric manifolds coincide. However, far from the Lagrangian points, the shapes of the apocentric manifolds allow to visualize the intricate chaotic dynamics known in dynamical systems’ terminology as the ‘homoclinic tangle’ (see Wiggins (1990)). Thus, the flux-tube and the apocentric manifolds are the same phase-space objects, but visualized differently in physical space. For a detailed comparison and generalizations of the various manifold theories see Harsoula et al. (2016), or the tutorial Efthymiopoulos (2010).

Since the manifold spirals arise in a system of reference which co-rotates with the bar, the manifold theory in its basic form predicts that the spiral arms should have the same pattern speed as the bar. This remark seems to come in conflict with observations both in our Galaxy (as reviewed e.g. in Bland-Hawthorn & Gerhard (2016); see also Antoja et al. (2014), Junqueira et al. (2015) and references therein) and in other galaxies (e.g. Vera-Villamizar et al. (2001); Boonyasait et al. (2005); Patsis et al. (2009); Meidt et al. (2009); Speights and Westpfahl (2012); Speights & Rooke (2016)). Regarding, now, galactic disc simulations, the leading paradigm over the years refers to simulations showing the co-existence of multiple pattern speeds (Sellwood & Sparke (1988); Little & Carlberg (1991); Rautiainen & Salo (1999); Quillen (2003); Minchev & Quillen (2006); Dubinski et al. (2009); Quillen et al. (2011); Minchev et al. (2012); Baba et al. (2013); Roca-Fabrega et al. (2013); Font et al. (2014); Baba (2015); but see also a noticeable exception in Roca-Fabrega et al. (2013)), possibly connected also to the phenomenon of nonlinear coupling of multiple disc modes (Tagger et al. (1987); Tagger & Athanassoula (1991); Sellwood & Wilkinson (1993); Masset & Tagger (1997)). On the other hand, it is well known that even isolated barred galaxies undergo substantial secular evolution (see Athanassoula (2013); Binney (2013); Kormendy (2013) in the tutorial volume Falcon-Barroso and Knapen (2013)). The tendency to transfer angular momentum outwards (e.g. towards the halo or across the disc, Tremaine & Weinberg (1984); Debattista & Sellwood (1998); Debattista & Sellwood (2000); Athanassoula (2002); Athanassoula & Misiriotis (2002); Athanassoula (2003); OŃeill & Dubinski (2003); Holley-Bockelmann et al. (2005); Berentzen et al. (2006); Martinez-Valpuesta et al. (2006)) leads the bar to slow down and grow in size at a rate which produces non-negligible change in dynamics at timescales comparable even to a few bar periods. This process becomes complex, and even partially reversed due to the growth of ‘pseudo-bulges’ or peanuts (Kormendy & Kennicutt (2004)), caused by dynamical instabilities such as chaos or the ‘buckling instability’ Combes & Sanders (1981); Combes et al. (1990); Pfenniger & Friedli (1991); Raha et al. (1991); Bureau & Athanassoula (1999); Martinez-Valpuesta & Shlosman (2004); Bureau & Athanassoula (2005); Debattista et al. (2006)). The reduction in size of the bar by the transfer of angular momentum under constant pattern speed is discussed in Weinberg & Katz (2007). Spiral activity acts as an additional factor of outwards transfer of angular momentum (Lynden-Bell & Kalnajs (1972)), while a radial re-distribution of matter can take place even under a nearly-preserved distribution of angular momentum (Hohl (1971); Sellwood & Binney (2002); Avila-Reese et al. (2005)). Radial migration is enhanced by the amplification of chaos due to the overlapping of resonances among the various patterns (Quillen (2003); Minchev & Quillen (2006); Quillen et al. (2011)).

As secular evolution affects both the bar’s morphology and pattern speed, as well as the extent to which other non-axisymmetric patterns are present in the disc, the shapes of the invariant manifolds found by momentarily ‘freezing’ the potential and pattern speed value also undergo important changes in time. Simulations by Lia Athanassoula have shown that, despite these changes, the stars out-flowing from the neighborhood of the Lagrangian points and develop orbits which, in general keep track of the change of the form of the invariant manifolds (Athanassoula (2012); see also Baba (2015); Lokas (2016)). As a rule, the material which populates the manifolds comes from orbital outflows originating from the interior of co-rotation, at the end of the bar (Contopoulos (1980)). As these outflows are adapted to the slowly-changing form of the manifolds, they are able to yield time-varying spiral or ring-like patterns. In principle, this non-rigidity of the shape and population density of the invariant manifolds leaves space to reconcile manifold spirals with time and space varying pattern speeds beyond the bar, although such a reconciliation awaits probe by specific simulations.

In the present paper we aim to investigate the application of the manifold theory of spirals in a typical model of secularly evolving bar. To this end, a series of isolated disc galaxy simulations were performed using the explicit Mesh-Adaptive N-Body method based on Approximate Inverses (MAIN, Kyziropoulos et al. (2016)), and the Parallel Self-Mesh Adaptive N-body technique based on Approximate Inverses (PMAIN, Kyziropoulos et al. (2017a); Kyziropoulos et al. (2017b)). The initial conditions for these simulations are described in Kyziropoulos et al. (2016). Whenever a bar is formed in a simulation, we find manifold spirals, as well as traces of secular evolution affecting the bar and other non-axisymmetric features in the disc. In the present paper we focus on one of these simulations (‘Q1’ in Kyziropoulos et al. (2016)), aiming to present these phenomena in considerable detail. This is a typical simulation starting with a Toomre exponential disc, a Dehnen-like ‘double-power-law’ halo, and a Sersic bulge, with parameters corresponding to a Milky-Way type galaxy (see Binney & Tremaine (2008)). The simulation undergoes phases of evolution known in literature since decades (see references above), namely: i) an initial instability leading to the formation of spirals, lasting for , ii) a bar instability, erasing all traces of the previous spirals and generating new ones related to the bar, iii) a ‘bucking instability’ which thickens the bar, and iv) a phase of slow evolution in which the amplitude of spirals is reduced. We focus our study on the epoch after the bar is formed (around Gyr), and up to . We measure several quantities which quantify the bar’s secular evolution. We then observe the corresponding evolution of the manifolds emanating from the neighborhood of and , as well as what happens, in general, to the disc as the bar-spiral structure evolves.

Our results lead to the following picture: as typical in such simulations, we find that spiral activity beyond the bar is characterized by recurrent episodes. ‘Incidents’ of growth of spiral and other non-axisymmetric disc features with substantial spectral power take place in a nearly damped-oscillatory manner. Damping, which results in these non-axisymmetric features gradually fading in amplitude, can be attributed to ‘disc heating’, i.e. the overall increase in time of the velocity dispersion in the outer parts of the disc (see Sellwood (2014) for a review). However, regarding what causes, in the first place, the repeated excitations of incidents of non-axisymmetric activity, besides recognizing in our simulation some well known mechanisms of enhancement of local perturbations (e.g. swing amplification, see also Baba (2015)), we find evidence of a mechanism not emphasized, to our knowledge, in literature. This is the fact that, every incident is accompanied by the creation of a population of particles moving in orbits either escaping or with apocenters far from the main part of the disc, whose distribution of directions of motion exhibits fluctuations with respect to perfect central symmetry. Hence, they cause the main part of the disc to recoil with respect to the spheroid (i.e. halo-bulge) center of mass. The recoil causes, in turn, an ‘off-centering effect’ resulting in a disturbance on the disc. Some studies have indicated the importance of such phenomena in the evolution of bar-halo systems, as well as the evolution of halo cusps (Sellwood (2003); McMillan & Dehnen (2005); Holley-Bockelmann et al. (2005); Debattista (2006); see also Hamilton et al. (2018) for a different application in the relaxation of globular clusters). Here we find that the ‘incidents’ of non-axisymmetric activity are correlated with the off-centering effect. The disc recoil sets the disc center of mass in relative orbit with respect to the spheroid (halo-bulge) center of mass. The relative orbit has small size (the ‘wobbling’ between disc and halo can extend as much as kpc, about of the bar size; off-centering effects of larger size are reported in McMillan & Dehnen (2005), with no apparent connection with the mechanism mentioned above). Yet, this is sufficient to induce a significant perturbation on the disc, with amplitude in a region including the bar. The disc response is found to have both and components. The latter might be connected to particles pushed to regions outside co-rotation, in which case they have to follow chaotic orbits guided by the (predominantly ) manifolds. Since at least one major recoil event is observed before the bar, the manifolds cannot be the drivers of the phenomenon. They can, however, help to maintain it, by supporting particles in chaotic orbits leading to fluctuations in the directions in which the particles move.

All together, we identify the role of the invariant manifolds as a skeleton of chaotic orbits in phase space, or, the dynamical avenues to be followed by particles whenever an incident of spiral or other non-axisymmetric is triggered in the disc. The morphology of observed patterns in the disc is always qualitatively similar to the ‘chaotic tangle’ formed by the manifolds, while direct comparison of the two patterns yields various levels of agreement, depending on the snapshot examined.

We finally study how all these phenomena evolve with the change in disc ‘temperature’ (velocity dispersion). The incidents of non-axisymmetric activity beyond the bar fade in a timescale Gyr. However, this happens in a purely stellar dynamical simulation. Invoking mechanisms able to dissipate kinetic energy in random motions (e.g. gas cooling, star formation etc.) should help to prolong these phenomena over a considerable fraction of the age of real disc galaxies.

The paper is structured as follows: section 2 presents the N-body simulation as well as methods of analysis. Section 3 explains the computation of the invariant manifolds and its comparison with the disc morphology. Section 4 deals with the disc’s secular evolution, and shows the comparison with the manifolds at various times. Section 5 summarizes our conclusions, along with some further comments on the role of manifolds in the modelling of barred-spiral galaxies.

2 Simulation and numerical computations

2.1 Initial conditions and N-body simulation

The simulation analyzed below corresponds to the experiment called ‘Q1’ in Kyziropoulos et al. (2016). In summary, we use particles to simulate:

i) an initially exponential disc of mass , exponential scale length kpc, vertical exponential scalelength kpc and Gaussian velocity distribution arising from a profile of Toomre’s Q-parameter rising in the center and tending to an asymptotic outward value with (see equation (13) of Kyziropoulos et al. (2016)).

ii) a Sersic-type spherical bulge embedded in the disc, mass , Sersic index , scale length kpc.

iii) A Dehnen-type double-power law spherical dark matter halo (see Binney & Tremaine (2008)), with density parameter , scale length kpc, and asymptotic inner and outer density power-law exponents , . The halo yields a mass at kpc, which rises to at kpc.

The disc, bulge and halo are represented by , and particles respectively. A ‘relaxation technique’ is used to match all three components in a unique system with virial ratio very close to (see Kyziropoulos et al. (2016) and Kyziropoulos et al. (2017a) for a description). Being prone to disc instabilities, the system exhibits vivid evolution dominated by the succession of spiral, bar, and other non-axisymmetric features.

The simulation was performed using the N-body code MAIN, which utilizes a Cartesian grid that solves Poisson equation by the use of a novel Symmetric Factored Approximate Sparse Inverse (SFASI)pre-conditioning technique, in conjunction with the multigrid method. Details on the features of the numerical technique are provided in Kyziropoulos et al. (2016). In Kyziropoulos et al. (2017a) and Kyziropoulos et al. (2017b) we presented two fully parallel versions of the N-body technique (for shared-memory and distributed-memory systems respectively). These versions implement fully adaptive mesh refinement while improving the algorithm by which boundary conditions are computed on the sides of the computational box. Using these improved techniques, we performed and particle versions of the numerical simulation Q1 (see Kyziropoulos et al. (2017a) and Kyziropoulos et al. (2017b)). The new simulations differ from the original one essentially only with respect to the exact moment when the bar instability is manifested. This is partly due to the better resolution of the profile of the central force, which affects the onset of the bar instability (Holley-Bockelmann et al. (2005)), and also to the fact that a smoother particle distribution implies lesser noise, and hence longer time for any instabilities to grow in the disc (e.g. by swing amplification, Toomre and Kalnajs (1991)). On the other hand, the main features of the bar and spirals (Fourier amplitudes, pattern speeds, secular evolution etc) appear similar in all runs.

Several computations require use of a smooth approximation of the potential in the disc plane . This is obtained, using bi-cubic interpolation, from the N-body potential evaluated in a grid , , where , kpc, .

2.2 Evolution of non-axisymmetric patterns

As a consequence of the chosen initial value of (=1), the disc exhibits a rapid growth of competing features. A spiral mode dominates up to about Gyr. After some transient phase, the bar starts growing rapidly between Gyr and Gyr. This growth is followed by ‘incidents’ of spiral activity (see below). Five such ‘incidents’ take place up to Gyr. The overall morphology of the galaxy changes rapidly in timescales as small as the bar period. The morphological evolution of the disc is shown in Fig. 1. As typical in such simulations, features like rings and spiral arms appear reccurently and with varying morphology after the bar formation. We focus on the incidents of spiral activity which take place between and Gyr. The bar performs slightly more than seven (counterclockwise) revolutions in this time interval. The spiral modes oscillate from conspicuous maxima (e.g. at , , Gyr) to very low minima (e.g. at , , , Gyr), albeit never vanishing completely. An inner ring-like structure surrounding the bar is also present, while, at particular snapshots (e.g. ) material appears to travel from one end of the bar to the other along the ring. This phenomenon is discussed in depth in several papers on the ‘flux-tube’ manifolds, e.g. Romero-Gomez et al. (2006); Athanassoula (2012).

The bar also thickens in time, forming a peanut, or ‘X-shaped’, pseudobulge). From an inspection of the vertical disc profiles and the evolution of the amplitude (see below), we infer that a buckling instability occurs around to Gyr. Manifold spirals, however, keep appearing well after Gyr. Thus, the manifold spirals do not appear to be halted by the buckling instability, as was reported for other simulations in Kwak et al. (2017) (see also Lokas (2016)).

Figure 1: Twenty snapshots of the face-on disc view in the interval . The bar rotates counter-clockwise. Spiral and ring-like structures appear reccurently and with varying amplitudes.
Figure 2: The non-axisymmetric excess density in nine snapshots in the interval . The color scale represents values in blue, in red and in between. The transition from a spiral to a bar mode dominance is shown. The first manifold spirals are developed immediately after the bar formation.
Figure 3: (a) Color map of the non-axisymmetric excess density at Gyr. A second incident of spiral activity accompanies the bar. (b) The Fourier amplitudes (cyan), (red), (green) and (blue). (c) The disc image by a logarithmic grayscale plot of the surface density . (d) The same image as in (c) after implementing the Sobel-Feldman edge detection algorithm (see text). Dimmer patterns of (c) are enhanced and clearly visible in (d).

We quantify the evolution of non-axisymmetric features using the disc surface density


where is the number of disc particles at time in each area element of a polar grid of 50 logarithmically equi-spaced radial bins from kpc to kpc, and 180 azimuthal bins from to . Next, we compute the Fourier transform


The relative amplitude and phase of the m-th Fourier mode are defined by


Truncating Eq. (2) in a finite number of harmonics () yields a smoothed representation of the surface density . The smoothed ‘non-axisymmetric excess density’ is defined as


The rapid growth of non-axisymmetric features in the disc leads to a dominant spiral mode lasting up to the time Gyr (Fig. 2). However, the onset of the bar instability erases all traces of previous spiral activity. The bar grows rapidly between Gyr and Gyr. A first conspicuous set of spiral arms beyond the bar appears around Gyr. The bar extends to kpc (as measured by the interval in in which the phase remains nearly constant), while co-rotation is initially at 6.5kpc. At Gyr, a second major incident of spiral activity takes place. Beyond this time, the and Fourier amplitudes enter a phase of considerably slower evolution.

Setting conventionally as the time at which the bar’s secular evolution begins, in the next section we focus on this time snapshot to describe the computations related to manifold spirals and the comparison between manifolds and observed N-body disc morphologies. As shown in Fig.3, the bar is very strong at this snapshot, with and at kpc. There is also important power in the and Fourier terms. The amplitude is larger than 0.1 across the bar’s whole extent and reaches in the central parts ( kpc). Such amplitudes are identifiable in several simulations (see, for example, Quillen et al. (2011); Minchev et al. (2012)). The question of what triggers the mode in the inner parts of the disc will be dealt with in section 4. On the other hand, in the snapshot of Fig. 3, substantial relative amplitudes are observed also at large distances ( kpc). This is connected to the lopesidedness of the observed spirals (and manifolds, see below), a fact which hints towards nonlinear interaction of the outer disc modes.

In Fig. 3, the density excess map (a) clearly shows the main spiral pattern, but not so clearly secondary features as the interarm features or the bifurcation of spirals seen in the lower left part of the direct image of the disc (panel (c)). In order to bring out these features, we produce an image of the disc (panel (d)) which enhances the edges of most ‘eye-recognizable’ patterns of image (c) using a pattern-recognition algorithm known as the ‘Sobel-Feldman edge detection’(see Gonzalez & Woods (1992)). This technique is similar to unsharp masking in that it processes the image through a smoothing filter (we used a Gausssian filter for image (c)). Instead of substracting the two images, Sobel-Feldman uses an operator to approximate the gradient of the smoothed image, and thus detect the edges of patterns, where the gradient becomes large. The Sobel-Feldman technique enhances the ability to recognize non-axisymmetric features of amplitude significantly lower than the one of the main bi-symmetric spirals. In subsequent sections we will use such ‘Sobel-Ferldman’ images of the disc in comparison with the patterns formed by the manifolds in the corresponding snapshots.

3 Manifold spirals

3.1 Definitions

At a fixed time , consider observers moving in a galactocentric frame rotating with angular velocity equal to the instantaneous angular velocity of the bar. Let be the total gravitational potential (arising from all components, e.g. bulge, disc-bar, halo) at the time . We consider first orbits in the rotating frame with ‘frozen’ in time potential . Let be the potential restricted in the disc plane. The orbits are given by Hamilton’s equations with Hamiltonian


where are cylindrical co-ordinates in the rotating frame, , , and are the radial and angular momentum (per unit mass) in the intertial frame, connected to the velocities in the rotating frame via , . The potential admits the Fourier decomposition:


Consider first a simplified bar-like model in which for all except and . The bar has a fixed orientation angle in the rotating frame. The equations of motion yield four Lagrangian equilibrium points (Binney & Tremaine (2008)). The unstable Lagrangian points are at the end of the bar, , . The stable Lagrangian points are at distance . The annulus is hereafter called the corotation zone. Introducing the full non-axisymmetric perturbation may induce an important deformation of the equipotential surfaces, influencing the number and position of unstable Lagrangian points (Tsoutsis et al. (2009), Athanassoula et al. (2009a), Athanassoula et al. (2009b), Kalapotharakos et al. (2010), Wu et al. (2016)). In the present simulation we observe only a small displacement of the points , with respect to the pure bar model.

The key property leading to the definition of the manifold spirals is now the following: one can show that there is a family of initial conditions for which the resulting orbits under the potential (6) tend asymptotically to the Lagrangian point when integrated backward in time, i.e., as . This leads to the formal definition of the unstable manifold of . Denote by one orbit with initial conditions , and also . The unstable manifold is the ensemble of all different initial conditions for which the distance between and tends to zero as :


Similar definitions hold for the unstable manifold of the Lagrangian point . Basic theorems of dynamics (Grobman (1959); Hartman (1960)) ensure that the sets , are non-empty. The dynamical role of the invariant manifolds and can be summarized as follows: collects the ensemble of trajectories tending arbitrarily close to in the backward sense of time. Thus, starting very close to , and integrating forward in time, these trajectories form an outflow away from . Projected in real space, this outflow yields the ‘flux-tube’ manifold (Romero-Gomez et al. (2006); Romero-Gomez et al. (2007)). An elementary analysis shows that the flux-tube manifold has two branches: one is directed outwards (outside co-rotation), and it takes the shape of a trailing spiral arm. The second is directed inwards (inside corotation), creating a ring-like structure around the bar. Thus, the flux-tube manifolds can give rise to several morphological structures, from rings to open spirals, depending on the model’s parameters (e.g. pattern spead, amplitude and asymmetry) as discussed in (Athanassoula et al. (2009a)).

In a similar way as for , we can define outflows of asymptotic trajectories emanating from a small but finite distance from . To this end, we invoke the family of short period (or ‘Lyapunov’) unstable periodic orbits which form small retrograde epicycles around , called hereafter the ‘PL1 orbits’ (Voglis et al. (2006a)) (denoted ). Similarly to , the unstable manifold of such an orbit, is defined as


The projection of the manifolds in the disc plane yields the ‘flux-tube’ manifolds . The geometric shapes and properties of the manifolds are similar to those of the manifolds . However, the periodic orbits PL1,2 and their manifolds form families which span a whole set of values of the Jacobi constant , while the unstable points and their manifolds represent only the values of the Jacobi constants . One has for any orbit of the PL1,2 families and their manifolds, implying that the manifolds characterize the motions of particles at all energies beyond the Jacobi energy at corotation.

The ‘flux-tube’ manifolds defined above represent a continuous swarm of trajectories with initial conditions in the set (3.1), or (3.1). This set reproduces the manifolds as geometric objects, but provides no information on the real distribution of matter along the manifolds, which depends on the distribution function of the N-body system. Despite the fact that a precise knowledge of the distribution function is hardly tractable, the method of apocentric manifolds (Voglis et al. (2006a)) exploits the fact that local density maxima along the manifolds are expected at points close to apsidal (i.e. pericentric or apocentric) positions of the orbits. This assumption is verified in N-body experiments (see Tsoutsis et al. (2008)), where it is shown that density maxima corresponding, in particular, to spiral patterns are associated with local apocenters of the orbits of the N-body particles. The apocentric manifolds are defined through the manifolds as follows:


A similar definition holds for the apocentric manifolds of the families PL1,2. For visualization purposes, the main benefit of the apocentric manifolds is that they allow to plot longer parts of the invariant manifolds, thus allowing to visualize in physical space the intricate chaotic dynamics induced by these manifolds, otherwise known as the ‘homoclinic’ or ‘chaotic tangle’ (see Efthymiopoulos (2010)). In the sequel we present computations based on the apocentric manifolds . In practical steps, we compute first an ‘apocentric surface of section’ for all trajectories of given Jacobi energy . Setting , , the equation


allows to compute the local radius at which a particle with trajectory corresponding to the Jacobi energy reaches a local apocentric passage with values of its angular variables equal to . Eq. (10) fixes the value of as a function of . Thus, for a randomly chosen trajectory undergoing epicyclic oscillations, the transition from one to the next apocentric passage can be viewed as a mapping


with the functions specified numerically via the integration of the trajectories. The orbits correspond to fixed points of the mapping (11), i.e., where the following condition is satisfied:


Eqs. (12) can be viewed as a algebraic system which can be solved numerically, using a root-finding technique, in order to compute the initial conditions of the corresponding PL1 or PL2 orbit. Numerical differentiation allows to compute also the elements of the monodromy matrix:


The matrix has two real eigenvalues satisfying . The linear eigenvectors associated with the absolutely largest eigenvalue define a slope, or ‘eigendirection’, of the linearized apocentric mapping in the neighborhood of the fixed point . To compute and visualize the apocentric invariant manifolds one takes many initial conditions equi-spaced along the unstable eigendirection and belonging to a linear segment of small total length on the apocentric surface of section . Integrating all these trajectories forward in time yields the manifolds . Taking only the iterates of the mapping (11) on the apocentric surface of section yields the apocentric manifolds . Note that the iterates of the mapping (11) are pairs of values . However, the constant energy condition (Eq. (10)) allows to compute the apocentric radius for any pair . Thus, any point on the apocentric surface of section corresponds to a triplet of values . The visualization of the apocentric manifold in physical space is obtained by plotting the manifold’s computed iterated points , . This is equivalent to computing the flux-tube manifolds and keeping only the points where the orbits have apocenter. On the other hand, the visualization of the apocentric manifolds in phase space (phase portrait) is obtained by plotting the manifold’s iterated points .

3.2 Manifold spirals in the N-body simulation

The computation of the apocentric manifolds in our N-body computation proceeds in the general way described above. Besides interpolating the N-body disc potential (see section 2), in order to compute the equations of motion, we first specify the value of the bar pattern speed as follows: computing the Fourier transform of the non-axisymmetric density excess (Eq. (2)), the mode is largely dominant for the bar. The angular displacement, for fixed radial distance , of the maxima of the mode at two successive snapshots separated by a time yield the pattern speed at the distance , namely , where is defined in Eq. (3) for . We use . We prefer this method over extending the Fourier transform (say of , ) in the time domain, since the latter approach assumes that patterns are characterized by constant frequencies in a relatively long time window, while the bar’s pattern speed appears to vary significantly over rather short time windows (see below). Figure 4 shows our calculation for the snaphsot Gyr. The curve exhibits an approximate ‘plateau’ at radii kpc. We estimate the bar pattern speed as the mean value of in the interval (the plateau is more clear in this range of radii). The vertical lines mark the positions of the Inner Linblad Resonance (ILR), co-rotation (CR), and outer Linblad resonance (OLR) estimated from the relations , , and , where is the epicyclic frequency derived from . The plateau is well formed outside the ILR radius, and yields a pattern speed km/sec/kpc. More comments on the structure of the curve are made in section 4.

Figure 4: The mode pattern speed for the snapshot Gyr. The bar pattern speed is estimated as the mean value of in the interval . The inner and outer dashed vertical lines mark the position of ILR and OLR. The solid vertical line marks the position of co-rotation (CR).

Having fixed we compute the unstable Lagrangian points and in the Hamiltonian (5). Since we use the full N-body potential, the Lagrangian points are not found at exactly symmetric positions with respect to the disc’s center, while their Jacobi energies have also a small difference . Similar differences hold for the stable Lagrangian points and . We finally use the information of the Lagrangian points to compute the apocentric manifold spirals emanating from the co-rotation zone. Depending on the choice of Jacobi energy, each orbit of the PL1 or PL2 families yields a different apocentric manifold. To end up with just one (representative) manifold calculation per family, we compute first the orbits PL1 and PL2 for various Jacobi energies, and keep, in each case, the orbit with Jacobi energy midway between the energies at and , i.e., , . These energies roughly correspond to the median of the Jacobi energy distribution for the N-body particles in the corotation zone. The apocentric unstable manifolds of the corresponding orbits PL1 and PL2 are then computed, taking an initial segment on the apocentric surface of section with N=10000 initial conditions distributed according to , , for the outer branch, and , for the inner branch of , with adjusted around a value so as to give comparable lengths for all computed manifolds, while , where is the slope of the unstable eigendirection of the the monodromy matrix at the fixed point PL1 (and similarly for the orbit PL2).

Figure 5: The apocentric surface of section as defined in Eq.(10) for the energy , at the snapshot Gyr. The orbits in this surface of section are strongly chaotic. The bottom panel is a magnification of the top panel in the region of angular momenta close to the value at corotation. The black thick points are the iterates of orbits with initial conditions taken at various locations within the chaotic domain. The red points show the the apocentric invariant manifolds . They emanate from the Lyapunov periodic orbits PL1 and PL2 shown as two blue fixed points (left and right respectively). The thick blue arrows indicate the corresponding unstable eigendirections. The closed green thick curves are limiting curves inside which the orbits are energetically forbidden. The marked features and are the ‘bridges’ and ‘gaps’ (see text).
Figure 6: Top: the apocentric manifolds (red) and (blue) superposed to the ‘Sobel-Feldman’ image of the disc at Gyr. Besides the spirals, the manifolds reproduce several secondary patterns in the disc beyond the bar, such as ‘gaps’, ‘bridges’, ‘rings’ and ‘bifurcations’ of secondary spiral arms (see text for details). Bottom: the apocentric manifolds together with the limiting curves of the apocentric surface of section (green), and with the bridges and gaps marked as in the corresponding surface of section of Fig. 5.

Figure 5 shows the form of the apocentric surface of section (plane ), as well as the apocentric manifolds of the PL1 and PL2 orbits, at the snapshot Gyr. The manifolds (red points) emanate as straight lines starting from the fixed points corresponding to the periodic orbits PL1 and PL2, but they soon develop a very complicated form, characterized by a number of thin lobes forming a ‘chaotic tangle’, which is the signature of homoclinic dynamics. Taking segments of initial conditions in the interior of the manifolds defined by the lobes marked , we iterate these chaotic orbits (black points), and find that the distribution of the iterates covers more and more uniformly the area inside the manifolds’ lobes.

All together, the phase portrait formed at the indicated level of energies (Jacobi constant ) has considerable chaos. Chaos can be partly due to numerical effects, (e.g. the ‘microchaos’ generated by noise in the potential computation), but it is also of dynamical origin, due to the strong bar and spiral amplitudes. In fact, while chaos fills a substantial part of the phase space at corotation already in models with strong bars (see, for example, Patsis et al. (1997)), the deformation of all phase space structures due to the spiral mode is evident in Fig. 5. This deformation affects the form of the closed limiting curves in the apocentric surface of section, inside which the motion is energetically forbidden. These curves are found by locating the pairs for which the effective potential satisfies the relation


where is the root of the equation . The above limiting curves are similar in shape to the usual curves of zero velocity, but provide the most stringent limits for the apocentric surface of section (see Tsoutsis et al. (2009).

In Fig. 5 the manifolds’ lobes from PL1 and PL2 encircle the corresponding limiting curves. However, despite the lack of energetic barrier, we observe that the uppermost manifold lobes do not fill the whole area around the limiting curves, but perform a number of oscillations leaving narrow gaps, denoted , in the surface of section. This behavior of the manifolds is dictated by basic rules of dynamics, which assert that the manifold lobes of the same or different periodic orbits cannot intersect each other. In this way, we observe that the lobes of the manifold emanating from PL2 approach the manifold emanating from PL1 in the area marked , but do not intersect it, leaving instead a gap . Similarly, the lobes of the manifold from PL1 approach the manifold from PL2 in the area , leaving a gap . These features have a specific morphological correspondence in physical space, as shown below. Let us note also that the level of Jacobi energy of the surface of section in Fig. 5 is close to the median of the entire distribution of the simulation’s particles in chaotic orbits, which is consistent with the general expectation that most particles supporting chaotic spirals should be distributed in the interval of energies between and (Patsis (2006); Tsoutsis et al. (2008))).

Figure 6 shows, now, the main result regarding the comparison between the disc morphology and the apocentric manifolds at the snapshot Gyr. The Sobel-Feldman image of the N-body disc shows a nearly bi-symmetric set of spirals, along with secondary ring and spiral features beyond the bar. The gaps and appear as narrow zones separating the manifolds at the bridges , . In these bridges, the manifold emanating from the region of approaches, in a nearly-tangent direction, the exterior side of the manifold emanating from the region of the Lagrangian point , and vice versa. The approach takes place via several oscillations of the manifolds in space, forming patterns recognized in the plot as bundles of preferential directions occupied by the manifolds. Such bundles mostly form spiral patterns, while breaks, or ‘bifurcations’, are also observed, splitting in two some of these bundles. Most notably, the main morphological features of the apocentric manifolds have counterparts in the Sobel-Feldman image of the real patterns formed by the disc particles.

As a final remark in this section, the manifolds and exhibit some lopesidedness, manifest also in the form of the limiting curves of the apocentric surface of section as projected in physical space. This effect implies that the outflow of particles in chaotic orbits, in the directions indicated by the manifolds, is not symmetric with respect to the disc center. Since the phase space outside corotation is open to escapes, particles escaping in chaotic orbits carry with them linear momentum in a distribution of orientations non-symmetric with respect to the center. The contribution of this to maintaining a disc-halo ‘off-centering effect’ is discussed in subsection 4.3.

4 Secular evolution

In this section we focus on how the various spiral and other non-axisymmetric features beyond the bar evolve along with the bar’s secular evolution. We also examine which mechanisms are responsible for the maintainance of appreciable levels of non-axisymmetric activity. In the end of the section we return to the comparison between manifolds and the observed patterns in the disc (as in Fig. 6), but for different snapshots covering the whole disc’s secular evolution.

4.1 Bar spin-down

A basic manifestation of secular evolution is the bar’s spin down, as shown in Fig. 7, which leads to a an outward slow displacement of the ILR, CR and OLR. Thus, kpc immediately after the bar formation and kpc at the end of the simulation. The bar’s spin down can be associated with several factors, as for example: i) particles abandoning the bar and carrying angular momentum outwards at incidences of growing spirals (see below), or ii) transfer of angular momentum to the halo (see Debattista & Sellwood (1998); Debattista & Sellwood (2000); Athanassoula (2002); Athanassoula & Misiriotis (2002); Athanassoula (2003)).

Figure 7: Top left: Evolution of the bar’s pattern speed. The cyan curve represents an exponential fitting , with km/sec/kpc, km/sec/kpc, Gyr, Gyr. The remaining panels show the evolution of the radii of the bar’s ILR (top right), CR (down left), and OLR (down right).

The bar retains high amplitudes throughout the simulation. Thus, whereas slowing down, the bar sustains appreciable levels of orbital chaos in the disc. Figure 8, top panel, shows the evolution of the distribution of the disc particles’ Jacobi energies, along with more phase portraits at energies, or snapshots, different from the one of Fig. 5. The distribution of Jacobi energies evolves by keeping its global maximum always close to the instantaneous CR value, hence it is shifted towards absolutely smaller energies as the bar slows down. The dispersion of the Jacobi energies is also reduced in time, but remains at a level of few times . The remaining panels in Fig. 8 show apocentric surfaces of section selected so as to cover the entire range of energies where particles are practically distributed, for Gyr (as in Fig. 5) and at a later snapshot Gyr. Similar phase portraits are found in all snapshots. Chaos is ubiquitous in all these portraits, indicating that the spiral and other structures beyond the bar are supported mostly by chaotic orbits.

Figure 8: Top left: Evolution of the distribution of the Jacobi energies of the disc particles at the snapshots Gyr (red), Gyr (green), and Gyr (blue). In comparison, the jacobi energies at CR are (in ), , and respectively. The top right and middle panels show the apocentric surfaces of section for the snapshot Gyr, in energies different from the one of Fig. 5, selected to cover the range of the corresponding distribution. The bottom panel shows two surfaces of section for the snapshot . Chaos is ubiquitous in all these plots.

4.2 Incidents of spiral and other non-axisymmetric activity

The variability of the bar, spiral and other non-axisymmetric patterns is reflected in time variations of the Fourier amplitude , as shown in Figure 9, top panel, for different radii. The group of curves for and kpc quantify the evolution of the bar’s strength. The bar reaches its maximum strength near Gyr, leading to a amplitude larger than unity. Thereafter, the bar amplitude stabilizes in the interval from to Gyr. After a small but abrupt drop at , associated with the buckling instability, the bar’s amplitude slowly decays through the remaining Gyr of the simulation, remaining however always of order of unity.

Figure 9: Top: evolution of the Fourier amplitude in the interval . The value of is shown at the distances kpc (dotted blue), kpc (solid line, blue), kpc (dotted green), and kpc (solid red). Bottom: evolution of the Fourier amplitude at kpc (dotted blue), kpc, (solid blue), and kpc (thick solid red).

The increase and subsequent slow decay of the bar mode permeates all radial annuli across the disc. The amplitude in general falls with between CR and OLR, but it tends to stabilize to a mean level in an annulus . The amplitude undergoes oscillations which are nearly in-phase for all radii within this annulus. At local maximum of these oscillations exceeds 0.2, up to Gyr, while the oscillations are milder afterwards. At all peaks of the mode in the interval we observe ‘incidents’ of non-axisymmetric activity in the outer parts of the disc. Taking the mean from the three lowermost curves in the top panel of Fig. 9, we estimate the times when such incidents reach maximum amplitude (to a precision related to the frequency of saving of the particles’ positions, i.e., every 25 Myr). We thus identify the following sequence of times:


By the above sequence, an approximate period Gyr can be deduced, but with large fluctuations Gyr. The bar has a pattern speed Km/sec/kpc (period Gyr) at the starting time and Km/sec/kpc ( Gyr) at the end of the sequence. Thus, the sequence (4.2) is in rough resonance with the bar’s rotation. However, the bar itself exhibits some extra phenomena. The bottom panel of Fig. 9 show important oscillations in the mode at disc radii smaller than the bar’s half-length. In these oscillations reaches a mean level , with peaks up to , and an approximate period Gyr in the interval from Gyr to Gyr. Oscillations of appear also in the outer disc, albeit not in phase with those of the inner disc.

Figure 10: Top: an ‘incident of inner origin’, whose corresponding maximum occurs around . Bottom: the Fourier amplitudes for (cyan), (red), (green), and (blue). The dotted vertical line marks the position of the bar’s CR. An disturbance propagates from CR outwards. Significant rise of the and modes appears inside the bar at the initial phase of the incident.

Through density excess maps, and the corresponding profiles of the Fourier modes , we can see how incidents of non-axisymmetric activity affect the disc’s morphology beyond the bar. We distinguish two types of incidents, exemplified in Figs.10 and 11 respectively, as follows:

-Incidents of inner origin: an incident of inner origin corresponds morphologically to a spiral wave emanating from the end of the bar, and travelling outwards, until it becomes detached from the bar. An example is given in Fig. 10, showing the density excess maps at nine snapshots around the time , which belongs to the sequence (4.2). The outwards ‘emission’ of the spiral wave from the bar becomes evident after Gyr. At earlier times, one observes the appearance of small leading extensions of the bar, so that the whole incident is reminiscent of swing amplification. The creation of a leading component can be associated with flow of material approaching the bar from the leading direction by moving along outer families of periodic orbits such as the 2:1 family (Patsis; private communication). The spirals start fading out after detaching from the bar.

Figure 11: Same as in Fig. 10, but for an incident of outer origin reaching its peak around . An disturbance generated in the outer parts of the disc is reflected at CR. Significant rise of the and modes appears again inside the bar.

- Incidents of outer origin: an incident of outer origin corresponds morphologically to a rise of the surface density excess in the outer parts of the disc, as shown in Fig. 11. The rise appears both in the and modes. The left panels in Fig. 11 indicate that the perturbation arises at radial distances beyond kpc, it travels inwards, and then it is reflected at corotation. It is stressed that these features are not connected with the present analysis of spirals based on manifolds. However, they may contribute to other effects observed in the disc as discussed below.

Figure 12 shows the density excess map of the disc for all times of the sequence (4.2). Incidents of inner origin appear at the snapshots , , , Gyr, while incidents of outer origin appear at , , , and Gyr. A classification as inner or outer is unclear for the incidents taking place at , , , and Gyr. Incidents of inner origin appear mostly in the interval , while incidents of outer origin are observed throughout the simulation. In fact, the disc becomes overall less responsive to non-axisymmetric perturbations as the time goes on, and a transition from stronger to weaker incidents is discernible around Gyr.

Figure 12: Non-axisymmetric density excess maps showing the peaks of the ‘incidents’ of non-axisymmetric activity for all snapshots of the sequence of Eq. (4.2).

4.3 Relative disc-halo orbit and its effect on the disc

At every inner incident of spiral activity, the particles’s orbits change as the potential develops rapid fluctuations. This is accompanied by changes in the particles’ energies and angular momenta. Figure 13 summarizes these effects. We consider moving time windows with Gyr, as well as an annulus , with kpc and kpc, equal to six disc exponential scale lengths. At a snapshot we isolate all particles whose orbits have crossed both limits of the annulus within the time window . We finally collect these particles for all time windows, and check their orbits. This process collects about 10% of the disc particles whose orbits behave as above at least in one time window. Figure 13, top, shows for a sample of 100 randomly chosen trajectories in this set. Analogous pictures are obtained for any other random choice. Note that the N-body code represents forces sufficiently smoothly to produce smooth particle orbits. The main feature of these orbits is their change of elongation (difference between apocenter and pericenter) in different time windows. Abrupt changes are observed for some orbits at particular times, connected to an abrupt change in the angular momentum (middle panel). These times correspond to the growth phase of major incidents of non-axisymmetric activity in the disc. Τhe effect is more notorious for the first incident of spiral growth a little before Gyr (before the bar), as well as the growth of the manifold spirals at and Gyr. In every incident the majority of affected orbits become more elongated after the incident than before (although both types or orbital changes are observed). Also, the majority of affected orbits exhibit increase of the angular momentum, which mostly affects an orbit by increasing its apocentric distance to a large value (more than 20 kpc for some particles in Fig. 13), while many particles move also to orbits with pericentric distance beyond co-rotation.

The relative influence of each incident in the orbits can be estimated by measuring the fraction , where is the total number of particles which, in at least one time window exhibited a crossing of the annulus as described above, and is the number of particles doing so in the time window . The bottom panel in Fig.13 shows the evolution of for four choices of , namely and kpc. Consecutive ‘incidents’ of non-axisymmetric activity cause an oscillatory behavior of these curves. For nearly all particles in the population the first event is associated with the growth of the bar (at Gyr), while for most particles (up to a percentage %) we have repeated events associated with the growth of the manifold spirals at , and Gyr. Smaller percentages appear at subsequent times, as the incidents of non-axisymmetric activity beyond the bar slowly decay in amplitude.

Figure 13: Top: evolution of the radius for one hundred randomly selected particles whose orbits cross the annulus in at least one time window (see text). The colored curves show three different types of behavior, namely, a gradual increase of the apocenter (red), an increase taking place mostly in isolated ‘incidents’ (blue), and an orbit going to a larger and then returning to a smaller apocenter (purple). Middle: evolution of the angular momentum per unit mass for the same particles. Bottom: the relative fractions (see text) when is taken equal to kpc (red), kpc (black), kpc (green), kpc (purple).
Figure 14: Percentage of disc particles found at a distance kpc, as a function of time. The arrows indicate times when an abrupt increase of escaping disc particles is observed. The dashed vertical lines correspond to the limits of time intervals ( to ) in which the relative disc-halo orbit is observed to change character, i.e., to switch from circumcentric to epicyclic, and vice versa (see text).

The orbital changes observed at major incidents of non-axisymmetric activity result in a increase in time of the percentage of disc particles beyond any sufficiently large radius . The increase is statistical, since it is due both to particles really escaping the system or ones moving in bounded, but elongated orbits, being pushed, however, outside co-rotation. We conventionally call the radius kpc as the ‘escape’ radius, and the particles found at at the time ‘escape bodies’. Figure 14 shows the evolution of the percentage of ‘escape bodies’ . The contribution of major incidents of non-axisymmetric activity is identified in this plot as an abrupt increase of . Three major events are observed to take place at the times , and Gyr, while smaller ones at the times and Gyr. Besides the first event at (before the bar), the remaining ones are all associated with incidents of inner origin. Outer incidents may also affect the particles’ orbits, but their influence on should be smaller, as the corresponding disc disturbances propagate inwards and then outwards (being reflected at CR, see above).

Figure 15: Top row: Time evolution of the coordinate (left), and coordinate (right) of the center of mass of the halo particles (blue dotted) or the bulge particles (cyan solid line). The thick green line shows the corresponding evolution of the and co-ordinates of the total center of mass of all particles included in the simulation box. Middle row: same as in the top row but for the centers of mass of the halo particles (blue-dotted) and the disc particles (red solid). Bottom row: Relative orbits of the halo (blue dotted) and disc (red solid) centers of mass with respect to the center of mass of the total system.

Assuming, now, that the particles connected with these events move in chaotic orbits, one expects important fluctuations in the number of particles escaping co-rotation as a function of the final direction in which the particles move. As a result, the main body of the disc recoils in a direction changing stochastically in time, the effect being larger at times connected with major incidents. This recoil sets the main body of the disc in relative orbit with respect to the spheroid components (i.e. bulge+halo). Figure 15 illustrates the phenomenon. We compute the evolution of the position of the center of mass , , , for the disc, bulge and halo components respectively. For the disc particles we use an iterative algorithm to specify the center of mass of only those particles characterized as ‘non-escaping’. Namely, we first compute the center of mass of all disc particles, and then consider the center of mass of only those disc particles satisfying , where is the cylindrical radial distance of a particle with respect to . Alternative ways to compute , taking the center of the computational box or the center of mass of the total system as origin, yield practically indistinguishable results. Finally, we compute the time evolution of the center of mass of all bodies in the simulation . In a perfect momentum-preserving scheme, the vector should remain constant in time (the net total momentum in the simulation initially vanishes). However, very slow smooth variations of the global center of mass are obtained as a result of the use of adaptive mesh-refinement as well as the fact that the accelerations of particles escaping the computational box are computed with multipole formulas instead of cloud-in-cell interpolation (see Kyziropoulos et al. (2017a) for details). These phenomena are a priori predictable, and they cause the global center of mass to undergo a smooth slow displacement with time-varying velocity not exceeding the value Km/sec in any of its three Cartesian components during the whole simulation (thick solid curves in top and middle panels of Fig. 15). On the other hand, the main effect to be observed in the same figure is that each of the three individual components , , undergo rapid stochastic variations around the global center of mass. We checked that such variations are important only in the projections of the three vectors on the disc plane, and they can be described as follows: The and components of the centers of mass of the halo and bulge undergo rapid stochastic variations around the smooth curves , , remaining, nevertheless always closely attached to each other. On the contrary, the and components of the vector of the disc center of mass exhibit rapid fluctuations which set the disc in relative orbit with the (nearly coinciding) centers of mass of the halo-bulge.

The evolution of the relative orbit of the halo-bulge, and disc with respect to the global center of mass is shown in the bottom panel of Fig. 15 ( components of and ). The separation between disc and halo-bulge increases near major incidents of non-axisymmetric activity, and the group of spheroidal components (halo-bulge) is set in relative motion with respect to the disc. The first major separation takes place near Gyr (first major growth of spirals before the bar). After the bar is formed, it is unclear whether the bar’s center of mass moves closer to the rest of the disc or to the spheroid (halo-bulge). We keep the conventional separation of the system into disc and spheroid (halo-bulge) components, hereafter referred to as the ‘relative disc-halo’ orbit . The disc-halo orbit is depicted in Fig. 16, starting from the time Gyr, immediately after the first incident of spiral activity accompanying bar formation. The four panels correspond to intervals in which the relative disc-halo orbit changes character, i.e., from ‘circumcentric’ (described around the disc center) the orbit switches to ‘epicyclic’ (the halo-bulge center of mass describes one or more small epicycles with respect to the disc center), and vice versa. These switches take place at the times Gyr, Gyr, Gyr, and Gyr. The first ‘epicyclic’ interval is rather short ( Gyr), while the second is longer ( Gyr) and within it the relative disc halo separation tends to zero, leading the whole system again close to a state of common center of mass of all components. Most notably, the times , , and differ by Gyr, from the times when abrupt changes in take place in Fig. 14. The time difference is of the order of one radial orbital period for the particles of Fig. 13. We regard these facts as evidence that the orbital switches observed in Fig. 16 are connected with the disc recoil effect.

Figure 16: The relative orbit of the spheroid (halo-bulge) center of mass with respect to a non-inertial observer centered at the disc center of mass (components vs. of the vector ) at four time intervals (top left), (top right), (bottom left), and (bottom right), with Gyr, Gyr, Gyr, Gyr.

The relative disc-halo orbit shown in Fig. 16 is quite small compared to the system’s scale. Fig. 17 shows the variations of from Gyr to Gyr. This interval can be split into sub-intervals representing circumcentric or epicyclic parts of the disc-halo orbit. In both cases, the distance presents an oscillatory behavior, leading to a minimum distance close to and maximum distance kpc. Despite the small size of the disc-halo orbit, the displacement of the halo-bulge component with respect to the disc causes a significant perturbation in the disc. Using a simple model (Appendix I) we find that the amplitude of the perturbation scales with the separation as at distances up to , where is the disc’s exponential scale length. With as small as kpc and kpc we find an induced perturbation (see Fig.22)). Figure 17 shows direct evidence of the correlation between variations of and of the relative amplitude in the disc. The correlation is strong in a domain (top panel), which contains the bar beyond its ILR, and connects maxima and minima of with maxima and minima of both in the circumcentric and epicyclic parts of the disc-halo orbit. Some change in the mean levels of both and is observed at every orbital switch. The correlation is weaker at larger radii (bottom panel). Thus, the displacement affects the disc’s inner part.

Figure 17: Top: evolution of the disc-halo displacement (red, values as in left vertical axis), compared to the mean from three values at and kpc (blue, values as in right vertical axis). The cyan vertical lines indicate the times to (same as in Fig. 16). Bottom: same as top but for the mean from three values at the radii and kpc.

The displacement causes also a direct effect whose size is of order , hence much smaller than the effect. However, the orbital response of disc particles to the perturbation affects mostly particles inside the bar. Due to chaos inside co-rotation (Fig. 8), even small potential fluctuations of any harmonics largely influence the particles’ orbits. In particular, particles coming closer to the manifolds can abandon co-rotation in the direction dictated by the manifolds. This leads to an indirect effect. The size of the effect can be estimated by the correlations between and the disc’s . Figure 18 shows that and tend to switch from anti-correlated to correlated when the disc-halo relative orbit switches from epicyclic to circumcentric. Uncertainties in this figure obstruct a clear conclusion. As additional evidence, Figs. 10 and 11 show a significant rise of the (and also ) component in the bar at the initial phase of the incidents of both inner and outer origin. A similar rise is observed in all major incidents.

Figure 18: Same as in 17, but for the mean (blue) at the indicated distances.

4.4 Manifold spirals and the evolution of non-axisymmetric patterns in the disc

It was argued above that the manifolds act as drivers of the orbital evolution for particles in chaotic orbits pushed to move away from co-rotation. This effect is pronounced at major incidents of non-axisymmetric activity, but the role of manifolds as drivers of the orbital evolution should hold independently of the growth or decay in amplitude of such incidents. We now compare the evolution of non-axisymmetric patterns in the disc beyond the bar with the evolution of the patterns created by the manifolds, both in minima and maxima of the non-axisymmetric activity. The comparison is in terms of the shapes of these patterns.

Figure 19: Left: Sobel-Feldman images of the disc in twelve snapshots corresponding to both maxima and minima of (see text). Right: the apocentric manifolds for the same snapshots (colors same as in Fig.6), superposed to the Sobel-Feldman images.

In Fig. 9, nearly in-phase oscillations of the amplitude occur in the outer parts of the disc. We focus on such oscillations in the time interval (major incidents of inner origin, connected with spiral patterns, do not appear beyond this time). To locate more precisely the times of occurence of such incidents, we choose a reference radius kpc which is between the bar’s CR and OLR throught the simulation. Six consecutive local maxima and six minima of occur at the times


Figure 19, left, shows the Sobel-Feldman images of the disc in the above times, alternating between a time of maximum and a time of minimum. In a timescale as short as one bar’s period one finds a strong variability of the patterns identified in the disc beyond the bar. Note that the Sobel-Feldman algorithm allows to recover patterns which are quite fuzzy in simple surface-density processed images of the disc. In particular, some form of spiral activity is identifiable in all snapshots of Fig. 19, although the morphology of these spirals (number of arms, pitch angle, radial extent, etc) rapidly change. Besides the main spirals, other secondary features appear and disappear in succession, as, for example, rings (e.g. at and Gyr) and bifurcations of more spirals (e.g. at Gyr). Figure 19, right, shows the same images on top of which the apocentric manifolds are superposed. An overall comparison shows that the manifolds and Sobel-Feldman deduced patterns exhibit a similar level of complexity (forming, in the case of the manifolds, a ‘chaotic tangle’ Wiggins (1990); see also Voglis et al. (2006a)). Also, they agree in shape and orientation to a large extent. In more detail, however, the level of agreement varies between snapshots. A (rather subjective) visual comparison, distinguishes snapshots of better (e.g. or Gyr) or worse agreement (e.g. Gyr) than average. The following are some basic remarks regarding this comparison:

- The manifolds’ homoclinic lobes account for the systematic appearance of features called gaps, bridges and bifurcations in section 3. Such features appear in nearly all panels of Fig. 19, and show a counterpart in the Sobel-Feldman disc images. We note that gaps and bridges can be recognized also in several plots of ‘flux-tube’ manifods (see Athanassoula (2012)), but they are better visualized using the apocentric manifolds.

- In nearly all patterns there are small phase differences ( between some manifold lobes and the corresponding spirals in the Sobel-Feldman image (the largest difference is found for Gyr). This might be a dynamical phenomenon (disc response has delay with respect to the manifolds), but it may also be an artefact of the visualization by the apocentric manifolds, since the disc’s local density maxima are close to, but do not necessarily coincide with the individual orbits’ apocenters (Tsoutsis et al. (2008)).

- For computational convenience, we just choose one value of the Jacobi constant for the manifold computation. The manifolds retain a rather robust pattern against changes of the Jacobi energy (Tsoutsis et al. (2008)), but small variations in shape are allowed. Most notably, besides the manifolds of the PL1 or PL2 orbits, the chaotic tangle is supported by the manifolds of all families of unstable periodic orbits in the corotation region. Superposed one on top of the other, these manifolds form amore robust spiral pattern called ‘manifold coalescence’ in Tsoutsis et al. (2008). The manifold coallesence is a better representation of chaotic dynamics in the co-rotation zone than the manifolds of any single family. Nevertheless, its computation for many snapshots is a voluminous work extending outside our present scope.

- Finally, despite the clear correlation between manifolds and Sobel-Feldman detected patterns, not all patterns in the disc need necessarily be linked to the manifolds. Structures such as outer rings or slowly rotating outer spirals may not be supported by chaotic flows, but instead, by the regular flows associated with standard density wave theory (Boonyasait et al. (2005); Struck (2015)).

4.5 Consistency with multiple pattern speeds

Figure 20: Profiles of the pattern speed as a function of the radius from the disc center at the same snapshots as in Fig. 19. The inner plateau at distances between kpc and kpc determines the bar’s pattern speed. The two vertical lines mark the position of CR and OLR. Snapshots close to a local maximum or minimum of the non-axisymmetric activity beyond the bar (parameterized by ; see text) are marked by ‘max’ and ‘min’ respectively.

The fact that multiple pattern speeds are detected in observations and simulations is considered as an argument against manifold-supported spirals (see, for example, Speights & Rooke (2016)). However, as argued above, measured pattern speeds can be affected both by the change in their morphology and by the transport of material along invariant manifolds (Athanassoula (2012)).

In order to quantify the evolution of pattern speeds, we produce profiles as in Fig. 4, but for different snapshots. Figure 20 shows these profiles for the same snapshots as in Fig. 19, corresponding to the consecutive maxima and minima of the amplitude. In all these snapshots, the inner plateau of nearly constant value marks the pattern speed of the bar. This plateau extends up to a distance equal to the bar’s CR. However, important fluctuations of the pattern speed beyond CR are observed at all snapshots after Gyr. From the fluctuating part of the profile between CR and the OLR, the following evolution pattern is identified: at the six snapshots of maximum of (‘max’ in Fig. 20), after a possible hump near corotation starts decreasing, but tends to stabilize again to a value as we approach the OLR. Measured values for this second ‘plateau’ are km/sec/kpc, which are a factor smaller than the (time decaying) value of . Such stabilization disappears at the minima of (‘min’ in Fig. 20), where constantly decreases with beyond CR. Most notably, the decrease leads to at a radius always close to the OLR. Decaying profiles of pattern speeds are reported in observations (e.g Speights and Westpfahl (2012); Speights & Rooke (2016)). It is of interest to check whether the criterion of where the decaying curve terminates can be exploited in real observations for the location of resonances. A systematic study of this topic is proposed.

On the other hand, Fig. 19 shows no appreciable difference in the levels of agreement between manifolds and disc patterns depending on whether we are at a maximum or minimum of . We conclude that the level of agreement is independent of the variability of pattern speeds beyond the bar. One may remark, in respect, that multiple patterns enhance chaos by the mechanism of resonance overlap (Chirikov (1979); Quillen (2003); Minchev & Quillen (2006)) and thus render more particles’ orbits ruled by manifold dynamics. On the other hand, whether or not the manifolds are populated by sufficiently many particles to dominate the global patterns in the disc depends on mechanisms able to inject new particles in chaotic orbits. Such mechanisms are distinct from the manifolds, as discussed in subsections 4.2 and 4.3.

4.6 Disc thermalization

As a final investigation, we examine the time evolution of the disc’s ‘temperature’, i.e., velocity dispersion profile, which is a crucial factor affecting both the disc’s responsiveness to perturbations as well as the particles’ orbits responsiveness to manifold dynamics (subsection 4.5). We quantify disc temperature in an annulus of width around the radius by the radial velocity dispersion , where is the radial velocity of the i-th particle in the annulus and . Similar results hold for the dispersion in the transverse and vertical velocity components.

Figure 21 shows the profile at four different times, namely , , and Gyr. At , the dependence of the radial velocity dispersion on can be approximated by the union of two exponential profiles, one for the inner part of the disc (in the interval ), with law , with Km/sec, kpc, kpc, and another for the outer part of the disc (), with law , with Km/sec, kpc.

Figure 21: Profile of the radial velocity dispersion at four different times, namely Gyr (thick blue with points), Gyr (thick green solid), Gyr (dotted red) and Gyr (solid red). The dashed vertical lines mark the position of CR, and the solid vertical lines mark the position of the OLR at Gyr (blue) and Gyr (red) respectively.

Three more curves in Fig. 21, corresponding to and Gyr, show how the profile versus evolves in time. At radii beyond the last position of the OLR ( kpc), the velocity dispersion increases in time by a factor between and in a time interval Gyr, and the corresponding exponential profile becomes less steep. We notice, however, that the profile becomes nearly horizontal, leading to an ‘isothermalization’ (constant velocity dispersion) in a domain of the disc roughly between kpc and kpc. As seen in Fig. 21, the inner limit of this domain nearly coincides with the innermost position of CR (at ), while the outer limit nearly coincides with the outermost position of the OLR (at Gyr) in the considered time interval. One remarks that, as they move outward, at least one of these resonances visits part of the interval , in which the isothermalization takes place. Physically, the domain is where chaotic motions completely dominate the dynamics. The corresponding particles belong to the well known ‘hot population’ (Sparke and Sellwood (1987)), whose orbits span the whole domain while recurrently entering also inside corotation. As a result, the strongly chaotic dynamics brings about an equalization of the velocity dispersion in the whole domain. On the other hand, we observe that the velocity dispersion profile remains practically invariant at distances smaller than the initial position of the ILR (at kpc). In fact, this resonance appears to act as a barrier for ‘heat’ (random kinetic energy) transfer across the disc.

The key remark, regarding the above evolution, is that although incidents of non-axisymmetric activity (in particular spirals) heat the outer parts of the disc, isothermalization in the domain scanned by the CR and OLR resonances implies that a part of the disc between CR and OLR becomes cooler with time. Physically, particles in the outer parts of the bar migrate outwards, carrying with them kinetic energy in the form of random motions. In Fig. 21, the initial and final profiles intersect at a radius kpc. The bar’s co-rotation at the time of the initial profile is at kpc, while it shifts to kpc at the final time. Estimating by the disc domain where manifolds rule the response of particles’ orbits to external perturbations, this domain has a width kpc at the initial time Gyr, which is limited to kpc at Gyr, while this domain is expected to shrink further to negligible widths as the bar’s co-rotation keeps moving outwards.

5 Discussion and conclusions

The origin and longevity of spiral structure in disc galaxies is still matter of intense debate. As new observational data and more accurate simulations come to surface, we gradually arrive at appreciating the potential complexity of the mechanisms which generate and maintain spiral structure. The co-existence and possible non-linear coupling of multiple patterns (see references in introduction) hints towards the chaotic nature of spiral arms, in particular when a strong bar is present. To the extent that spirals are bar driven, they are affected in a non-trivial way by the bar’s secular evolution. Within this context, the main contribution of the manifold theory is to describe what should be the expected form of the bar-driven spiral mode when the disc’s region between the bar’s Corotation and Outer Lindblad Resonance is largely chaotic. It should be stressed that the manifold theory poses no requirement that chaos originates exclusively from the bar. Nonlinear interaction of the bar mode with additional patterns beyond co-rotation actually enhances chaos (see references in text). Furthermore, it is a basic rule of dynamics that the unstable manifolds of one periodic orbit can intersect neither themselves nor the unstable manifolds of any other periodic orbit of equal Jacobi energy. Thus, all these manifolds have to ‘coalesce’ in nearly parallel directions, thus enhancing chaotic spirals.

The manifolds emanating from the region of the bar’s and points provide the simplest representation of these chaotic outflows. In the present paper we give evidence that these manifolds provide a dynamical skeleton in phase space, or, the dynamical avenues to be followed by new particles injected in the domain between CR and OLR at reccurent ‘incidents’ of non-axisymmetric activity. We gave a characterization of such incidents as of i) inner, or ii) outer origin, depending on whether the spectral analysis shows a wave originating inside CR and propagating outwards (in (i)), or originating outside CR, moving initially inwards, then being reflected at CR, and then moving outwards (in (ii)). Manifold spirals are connected with incidents of inner origin. Pattern detection algorithms such as Sobel-Feldman allow to detect the co-existence of various patterns, both at maxima and minima of the amplitude beyond the bar. We interpret the importance of manifolds as follows: whatever causes these incidents, at every incident the particles’ orbits (and in particular chaotic ones) are perturbed. Then, new particles injected in the CR zone tend to follow and populate the manifolds to a large extent, according to general rules of dynamics.

We argued that the continuous change of the form of the manifolds, as well as the motion, along the manifolds, of the matter which populates them, allows to reconcile manifold spirals with the multiplicity (and variability) of pattern speeds beyond the bar. A basic behavior is found to govern the evolution of the radial profile of the pattern frequency beyond the bar. At maxima of the non-axisymmetric activity tends to form a second plateau, indicating a second pattern speed outside CR, well distinct from . At minima, the second plateau disappears, giving its place, to a shear-indicating decaying profile of the curve . Remarkably, in the latter case, terminates at a value always near . Five full cycles of this behavior are seen in our simulation, in a period of Gyr (see Fig. 20), leading to an approximate period of Gyr. This is in rough resonance with the bar’s period, but with large uncertainties in the numbers. On the other hand, since the comparison between manifolds and Sobel-Feldman-recognized patterns in the disc (Fig. 19) shows no appreciable differences between snapshots of minima and maxima of the non-axisymmetric activity, we argued that the degree to which manifolds are able to dictate the dynamics outside CR seems to be rather independent of the strength of any additional pattern in the disc.

In addition to well known mechanisms, we discussed a new mechanism which seems to play a key role in triggering new incidents of non-axisymmetric activity: this is the relative halo-disc orbit which results from the off-centering effect: whenever particles are ejected away from the disc through manifold spirals, the disc rebounces, leading to a relative orbit between the disc and the remaining spheroid (halo+bulge) component of the galaxy. We argue that even moderate off-centerings, of order , where is the off-center displacement and the disc’s exponential scale length, are able to induce appreciable tides on the disc (of order ). The corelation between and the disc response is clearly seen in the disc’s spectral analysis. The way the tide affects the response beyond the bar appears to be a complex nonlinear phenomenon for which we have no precise analysis at present. Yet, the fact that one phenomenon influences the other is clearly seen in our simulation, as detailed in subsections 4.2 to 4.4.

As a final remark, we discussed the ‘thermal’ evolution of the disc, i.e., the way in which the radial profile of the velocity dispersion appears to be influenced in time due to non-axisymmetric activity. Owing, again, to the large degree of chaos between CR and OLR, we argued that the gradual outward shift of both the CR and OLR radii as the bar slows down, causes part of the disc, beyond the end of the bar, to acquire a nearly constant radial velocity dispersion. In fact, this ‘isothermalization’ causes a part of the disc, starting from inside the bar and ending at a point midway between CR and OLR to cool down. We interpret this effect as a hint that chaotic populations of particles gradually migrate outwards, carrying with them kinetic energy in the form of randomly oriented motions. Regarding, however, the responsiveness of the disc to manifold dynamics, we provide a heuristic argument showing that good levels of responsiveness should be limited in a domain which shrinks in time, as the CR radius moves outwards, while the radius beyond which the disc gets hotter in time is nearly fixed.


We thank the reviewer, Prof. J. Binney, for his detailed report with many suggestions for improvement. We thank Prof. E. Athanassoula for a very careful reading of the manuscript, with many critical comments, as well as for enlightening discussions on manifold dynamics. We also thank Prof. P. Patsis for useful discussions. The authors acknowledge the Greek Research and Technology Network (GRNET) for the provision of the National HPC facility ARIS under Project PR004033-ScaleSciComp II. R.I. Paez was funded by the Project 200/854 of the Research Committee of the Academy of Athens and the ERC Project 677793 StableChaoticPlanetM, in non overlapping periods. The potential interpolation algorithm and manifold computations were developed as part of K. Zouloumi’s MSc thesis (Department of Physics, University of Athens, 2017).


  • Antoja et al. (2014) Antoja, T, Helmi, A., Dehnen, W., et al.: 2014, A&A, 563, 60
  • Athanassoula (2002) Athanassoula, E.: 2002, ApJ, 569, L83
  • Athanassoula (2003) Athanassoula, E.: 2003, MNRAS, 341, 1179
  • Athanassoula (2007) Athanassoula, E.: 2007, MNRAS, 377, 1569
  • Athanassoula & Misiriotis (2002) Athanassoula, E., and Misiriotis, A.: 2002, MNRAS, 330, 35
  • Athanassoula et al. (2009a) Athanassoula, E., Romero-Gómez, M., Masdemont, J. J.: 2009a, MNRAS, 394, 67
  • Athanassoula et al. (2009b) Athanassoula, E., Romero-Gómez, M., Bosma, A., Masdemont, J.J.: 2009b, MNRAS, 400, 1706
  • Athanassoula et al. (2010) Athanassoula, E., Romero-Gómez, M., Bosma, A., Masdemont, J.J.: 2010, MNRAS, 407, 1433
  • Athanassoula (2012) Athanassoula, E.: 2012, MNRAS, 426, L46
  • Athanassoula (2013) Athanassoula, E.: 2013, in ‘Secular Evolution of Galaxies’, XXIII Canary Islands Winter School of Astrophysics. J Falcon-Barroso and J. H. Knapen (eds), Cambridge University Press, pp.305-352.
  • Avila-Reese et al. (2005) Avila-Reese, V., Carrillo, A., Valenzuela, O., and Klypin, A.: 2005, MNRAS, 361, 997
  • Baba et al. (2013) Baba, J., Saitoh, T.R., Wada, K.: 2013, ApJ, 763, 46
  • Baba (2015) Baba, J.: 2015, MNRAS, 454, 2954
  • Berentzen et al. (2006) Berentzen, I., Shlosman, I., and Jogee, S.: 2006: ApJ, 637, 582
  • Binney & Tremaine (2008) Binney J., and Tremaine S.: 2008, Galactic Dynamics, second edn. Princeton University Press
  • Binney (2013) Binney, J.: 2013, in ‘Secular Evolution of Galaxies’, XXIII Canary Islands Winter School of Astrophysics. J Falcon-Barroso and J. H. Knapen (eds), Cambridge University Press, pp.259-304.
  • Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J., and Gerhard, O.: 2016, Annu. Rev. A&A Astrophys., 54, 529
  • Boonyasait et al. (2005) Boonyasait, V., Patsis, P.A., and Gottesman, S.T: 2005, NYASA, 1045, 203
  • Bureau & Athanassoula (1999) Bureau, M., and Athanassoula, E.: 1999, ApJ, 522, 686
  • Bureau & Athanassoula (2005) Bureau, M., and Athanassoula, E. 2005, ApJ, 626, 159
  • Chirikov (1979) Chirikov, B.V.: 1979, Phys. Rep., 52, 263
  • Combes & Sanders (1981) Combes, F., and Sanders, R.H.: 1981, A&A, 96, 164
  • Combes et al. (1990) Combes, F., Debbasch, F., Friedli, D., and Pfenniger, D.: 1990, A&A, 233, 82
  • Contopoulos (1980) Contopoulos,G.: 1980, A&A, 81, 198
  • Contopoulos (2002) Contopoulos,G.: 2002, ‘Order and Chaos in Dynamical Astronomy’, Springer
  • Danby (1965) Danby, J.M.A., 1965: AJ, 70, 501
  • Debattista & Sellwood (1998) Debattista, V.P., and Sellwood, J.A.: 1998, ApJ, 493, L5
  • Debattista & Sellwood (2000) Debattista, V.P., and Sellwood, J.A.: 2000, ApJ, 543, 704
  • Debattista et al. (2006) Debattista, V. P., Mayer, L., Carollo, C. M., Moore, B., Wadsley, J., and Quinn, T.: 2006, ApJ, 645, 209
  • Debattista (2006) Debattista, V.P.: 2006, ASPC, 352, 161
  • Dobbs and Baba (2014) Dobbs, C. and Baba, J.: 2014, Pub. Astron. Soc.