Escape dynamics and fractal basin boundaries in Seyfert galaxies
Abstract
The escape dynamics in a simple analytical gravitational model which describes the motion of stars in a Seyfert galaxy is investigated in detail. We conduct a thorough numerical analysis distinguishing between regular and chaotic orbits as well as between trapped and escaping orbits, considering only unbounded motion for several energy levels. In order to distinguish safely and with certainty between ordered and chaotic motion, we apply the Smaller ALingment Index (SALI) method. It is of particular interest to locate the escape basins through the openings around the collinear Lagrangian points and and relate them with the corresponding spatial distribution of the escape times of the orbits. Our exploration takes place both in the physical and in the phase space in order to elucidate the escape process as well as the overall orbital properties of the galactic system. Our numerical analysis reveals the strong dependence of the properties of the considered escape basins with the total orbital energy, with a remarkable presence of fractal basin boundaries along all the escape regimes. It was also observed, that for energy levels close to the critical escape energy the escape rates of orbits are large, while for much higher values of energy most of the orbits have low escape periods or they escape immediately to infinity. We also present evidence obtained through numerical simulations that our model can describe the formation and the evolution of the observed spiral structure in Seyfert galaxies. We hope our outcomes to be useful for a further understanding of the escape mechanism in galaxies with active nuclei.
Keywords:
Hamiltonian systems; galaxies; numerical simulations; escapes; fractals∎
1 Introduction
Over the last decades a huge amount of research work has been devoted on the subject of escaping particles from open dynamical systems. Especially the issue of escape in Hamiltonian systems is a classical problem in nonlinear dynamics (e.g., [C90, , CK92, – CHLG12, , STN02, ]). It is well known, that several types of Hamiltonian systems have a finite energy of escape. For values of energy lower than the escape energy the equipotential surfaces of the systems are close which means that orbits are bound and therefore escape is impossible. For energy levels above the escape energy on the other hand, the equipotential surfaces open and exit channels emerge through which the particles can escape to infinity. The literature is replete with studies of such “open" Hamiltonian systems (e.g., [BBS09, , KSCD99, , NH01, , STN02, , SCK95, – SKCD96, , Z13, , Z14, ]). At this point we should emphasize, that all the abovementioned references on escapes in Hamiltonian system are exemplary rather than exhaustive, taking into account that a vast quantity of related literature exists.
Nevertheless, the issue of escaping orbits in Hamiltonian systems is by far less explored than the closely related problem of chaotic scattering. In this situation, a test particle coming from infinity approaches and then scatters off a complex potential. This phenomenon is well investigated as well interpreted from the viewpoint of chaos theory (e.g., [BGOB88, – BOG89, , BM92, , CDO90, , DGO90, – EJ86, , GR89, , H88, , JRS92, , J87, – JT91, , LMG00, – LJ99, , ML02, , RJ94, , SASL06, – SS10, ]). Chaotic scattering has also been applied in the astrophysical context of many aspects such as scattering off black holes (e.g., [A97, , dML99, ]) and threebody stellar systems (e.g., [BTS96, , BST98, , H83, , HB83, ]). The related invariant manifolds of the chaotic saddle is directly associated with the chaotic dynamical behavior (e.g., [O02, ]). In particular, the chaotic saddle is defined as the intersection of its stable and unstable manifolds [S14, ], while hyperbolic and nonhyperbolic chaotic saddles may occur in dynamical systems (e.g., [LGB93, ]). More details on the issue of chaotic scattering and escape from chaotic systems can be found to the recent reviews [SS13, ] and [APT13, ], respectively.
In open Hamiltonian systems an issue of paramount importance is the determination of the basins of escape, similar to basins of attraction in dissipative systems or even the NewtonRaphson fractal structures. An escape basin is defined as a local set of initial conditions of orbits for which the test particles escape through a certain exit in the equipotential surface for energies above the escape value. Basins of escape have been studied in many earlier papers (e.g., [BGOB88, , C02, , KY91, , PCOG96, ]). The reader can find more details regarding basins of escape in [C02, ], while the review [Z14ip, ] provides information about the escape properties of orbits in a multichannel dynamical system of a twodimensional perturbed harmonic oscillator. The boundaries of an escape basins may be fractal (e.g., [AVS09, , BGOB88, ]) or even respect the more restrictive Wada property (e.g., [AVS01, ]), in the case where three or more escape channels coexist in the equipotential surface.
One of the most characteristic Hamiltonian systems of two degrees of freedom with escape channels is undoubtedly the wellknown HénonHeiles system [HH64, ]. A huge load of research on the escape properties of this system has been conducted over the years (e.g., [AVS01, – AVS03, , BBS08, , BBS09, ]). Escaping orbits in the classical RestrictedThreeBodyProblem (RTBP) is another typical example (e.g., [N04, , N05, , TdA14, ]). Furthermore, escaping and trapped motion of stars in stellar systems is an another issue of great importance. In a recent article [Z12a, ], we explored the nature of orbits of stars in a galactictype potential, which can be considered to describe local motion in the meridional plane near the central parts of an axially symmetric galaxy. It was observed, that apart from the trapped orbits there are two types of escaping orbits, those which escape fast and those which need to spend vast time intervals inside the equipotential surface before they find the exit and eventually escape. Furthermore, the chaotic dynamics within a star cluster embedded in the tidal field of a galaxy was explored in [EJS08, ]. In particular, by scanning thoroughly the phase space and obtaining the basins of escape with the respective escape times it was revealed, that the higher escape times correspond to initial conditions of orbits near the fractal basin boundaries.
The aim of this work is to numerically investigate the escape dynamics in a Seyfert galaxy with a rotating active central nucleus. The recent paper of Ernst & Peters (2014) [EP14, ] is used as a guide since we follow the same numerical methods and computational approach. The present paper however, contains several new aspects with respect to [EP14, ] that can be considered as further developments in the field of open galactic systems. In particular: (i) we proceed one step further by using the SALI method in order to distinguish between regular and chaotic trapped motion, (ii) we display the evolution of the percentages of all types of orbits in both the physical and phase space as a function of the value of the energy integral, (iii) we conduct an overview analysis thus revealing how the orbital structure of the dynamical system is influenced by the total energy when we use a continuous spectrum of energy values rather than some few discrete levels and (iv) we relate the escape times of orbits with the energy as well as with the width of the symmetrical escape channels.
The article is organized as follows: In Section 2 we present in detail the structure and the properties of our Seyfert galaxy model. All the computational methods we used in order to determine the character of orbits are described in Section 3. In the following Section, we conduct a thorough numerical investigation revealing the overall orbital structure (bounded regions and basins of escape) of the galaxy and how it is affected by the total orbital energy. In Section 5 we show that our dynamical model can simulate, in away, the creation as well as the evolution of spiral arms of a real Seyfert galaxy. Our paper ends with Section 6, where the discussion and the conclusions of this research are given.
2 Theory of the galaxy model
For the description of the Seyfert galaxy we use the following massmodel potential
(1) 
where is the gravitational constant, is the total mass of the galaxy, is a parameter leading to deviation from axial and spherical symmetry, while is the scale length of the galaxy which acts also as a softening parameter. The generalized anharmonic Plummer potential (1) has been successfully applied several times in the past in order to realistically model non axially symmetric galaxies (see e.g., [CP07, , MAB95, , PC06, , Z12b, ]). At this point we would like to stress out that the singlecomponent potential of Eq. (1) models the properties of the Seyfert galaxy as a whole rather than an isolate component (e.g., the galactic bar).
We consider the case where the dense and compact active nucleus of the Seyfert galaxy rotates clockwise at a constant angular velocity . Therefore, the effective potential is given by
(2) 
which is in fact a cut through the equatorial plane of the threedimensional (3D) potential.
In our study, we use a system of galactic units, where the unit of length is 1 kpc, the unit of mass is and the unit of time is yr (approximately 100 Myr). The velocity unit is 10 km/s, while is equal to unity . The energy unit (per unit mass) is 100 kms. In these units, the values of the involved parameters are: (corresponding to 1.1625 ), , , and they remain constant throughout the numerical analysis. The values of and were chosen having in mind the Seyfert galaxy NGC 1566 (see at the end of section 5 for more details).
This galactic model has five equilibria called Lagrangian points [S67, ] at which
(3) 
The isolines contours of constant effective potential as well as the position of the five Lagrangian points are shown in Fig. 1a. Three of them, , , and , known as the collinear points, are located in the axis. The central stationary point at is a local minimum of . At the other four Lagrangian points it is possible for the test particle to move in a circular orbit, while appearing to be stationary in the rotating frame. For this circular orbit, the centrifugal and the gravitational forces precisely balance. The stationary points and at are saddle points, where is called Lagrangian radius. Let located at , while be at . The points and on the other hand, which are located at are local maxima of the effective potential, enclosed by the bananashaped isolines. The annulus bounded by the circles through , and , is known as the “region of coroation" (see also [BT08, ]).
The equations of motion in the rotating frame are described in vectorial form by
(4) 
where is the position vector, is the constant rotation velocity vector around the vertical axis, while the term represents the Coriolis force. Decomposing Eq. (4) into rectangular Cartesian coordinates , we obtain
(5) 
where the dot indicates derivative with respect to the time.
In the same vein, the equations describing the evolution of a deviation vector which joins the corresponding phase space points of two initially nearby orbits, needed for the calculation of standard chaos indicators (the SALI in our case) are given by the following variational equations
(6) 
Consequently, the corresponding Hamiltonian (known also as the Jacobian) to the effective potential given in Eq. (2) reads
(7) 
where and are the momenta per unit mass, conjugate to and , respectively, while is the numerical value of the Jacobian, which is conserved since it is an isolating integral of motion. Thus, an orbit with a given value for it’s energy integral is restricted in its motion to regions in which , while all other regions are forbidden to the star. The numerical value of the effective potential at the two Lagrangian points and , that is and , respectively yields to a critical Jacobi constant , which can be used to define a dimensionless energy parameter as
(8) 
where is some other value of the Jacobian. The dimensionless energy parameter makes the reference to energy levels more convenient. For the equipotential surface encloses the critical area^{1}^{1}1In three dimensions the last closed equipotential surface of the effective potential passing through the Lagrangian points and encloses a critical volume., while for a Jacobian value , or in other words when , the equipotential surface is open and consequently stars can escape from the galaxy. In Fig. 1b we present the contour of when . We observe, the two openings (exit channels) through which the particles can leak out. In fact, we may say that these two exits act as hoses connecting the interior region (green color) of the galaxy where with the “outside world" of the exterior region (amber color). Exit channel 1 at the negative direction corresponds to escape through Lagrangian point , while exit channel 2 at the positive direction indicates escape through . The forbidden regions of motion within the bananashaped isolines around and points are shown in Fig. 1b with gray. It is evident, that for the two forbidden regions merge together thus escape is impossible since the interior region cannon communicate with the exterior area.
3 Computational methods
In order to explore the escape dynamics of stars in the Seyfert galaxy model, we need to define samples of initial conditions of orbits whose properties (bounded or escaping motion) will be identified. The best method for this purpose, would have been to choose the sets of initial conditions of the orbits from a distribution function of the system. This however, is not available, so we define for each value of the total orbital energy (all tested energy levels are above the escape energy), dense grids of initial conditions regularly distributed in the area allowed by the value of the energy. Our investigation takes place both in the physical and in the phase space for a better understanding of the escape process. In both cases, the step separation of the initial conditions along the axes (or in other words, the density of the grids) is controlled in such a way that always there are about orbits (a maximum grid of 318 318 equally spaced initial conditions of orbits), depending on the value of the Jacobian integral. Furthermore, the grids of initial conditions of orbits whose properties will be examined are defined as follows: For the physical space we consider orbits with initial conditions with , while the initial value of is always obtained from the energy integral (7) as . Similarly, for the phase space we consider orbits with initial conditions with , while this time the initial value of is obtained from the Jacobi integral (7).
The equations of motion as well as the variational equations for the initial conditions of all orbits were integrated using a double precision BulirschStoer FORTRAN 77
algorithm (e.g., [PTVF92, ]) with a small time step of order of , which is sufficient enough for the desired accuracy of our computations. Here we should emphasize, that our previous numerical experience suggests that the BulirschStoer integrator is both faster and more accurate than a double precision RungeKuttaFehlberg algorithm of order 7 with CashKarp coefficients. Throughout all our computations, the Jacobian energy integral (Eq. (7)) was conserved better than one part in , although for most orbits it was better than one part in .
In dynamical systems with escapes an issue of paramount importance is the determination of the position as well as the time at which an orbit escapes. For all values of energy smaller than the critical value (escape energy), the equipotential surface is closed and therefore all orbits are bounded. On the other hand, when the equipotential surface is open and extend to infinity. The value of the energy itself however, does not furnish a sufficient condition for escape. An open equipotential curve consists of two branches forming channels through which an orbit can escape to infinity (see Fig. 1b). It was proved [C79, ] that in Hamiltonian systems there is a highly unstable periodic orbit at every opening close to the line of maximum potential which is called a Lyapunov orbit. In fact, there is a family of unstable Lyapunov orbits around each collinear Lagrangian point [M58, ]. Such an orbit reaches the Zero Velocity Curve^{2}^{2}2The boundaries of the accessible regions of motion are the called Zero Velocity Curves, since they are the locus in the space where kinetic energy vanishes. (ZVC), on both sides of the opening and returns along the same path thus, connecting two opposite branches of the ZVC. Lyapunov orbits are very important for the escapes from the system, since if an orbit intersects any one of these orbits with velocity pointing outwards moves always outwards and eventually escapes from the system without any further intersections with the surface of section (e.g., [C90, ]. When the Lagrangian points exist precisely but when an unstable Lyapunov periodic orbit is located close to each of these two points (e.g., [H69, ]). The two unstable Lyapunov orbits around and are shown in Fig. 1b in red.
These two Lagrangian points are saddle points of the effective potential, so when , a star must pass close enough to one of these points in order to escape. Thus, in our galactic system the escape criterion is purely geometric. In particular, escapers are defined to be those stars moving in orbits beyond the Lagrangian radius thus passing one of the two Lagrangian points ( or ) with velocity pointing outwards. Here we must emphasize that orbits with initial conditions outside the interior region, or in other words outside or , does not necessarily escape immediately from the galaxy. In Fig. 2 we present a characteristic example of such an orbit with initial conditions: , , when . We observe that even though this orbit is initiated outside but very close to it does not escape right away from the system. On the other hand, it enters the interior galactic region and only after 32 time units of chaotic motion it eventually escapes from the opposite channel. Thus it is evident, that the initial position itself does not furnish a sufficient condition for escape, since the escape criterion is a combination of the coordinates and the velocity of stars.
The physical and the phase space are divided into the escaping and nonescaping (trapped) space. Usually, the vast majority of the trapped space is occupied by initial conditions of regular orbits forming stability islands where a third integral is present. In many systems however, trapped chaotic orbits have also been observed. Therefore, we decided to distinguish between regular and chaotic trapped motion. Over the years, several chaos indicators have been developed in order to determine the character of orbits. In our case, we chose to use the Smaller ALingment Index (SALI) method. The SALI [S01, ] has been proved a very fast, reliable and effective tool, which is defined as
(9) 
where and are the alignments indices, while and , are two deviation vectors which initially point in two random directions. For distinguishing between ordered and chaotic motion, all we have to do is to compute the SALI along time interval of numerical integration. In particular, we track simultaneously the timeevolution of the main orbit itself as well as the two deviation vectors and in order to compute the SALI.
The timeevolution of SALI strongly depends on the nature of the computed orbit since when an orbit is regular the SALI exhibits small fluctuations around non zero values, while on the other hand, in the case of chaotic orbits the SALI after a small transient period it tends exponentially to zero approaching the limit of the accuracy of the computer . Therefore, the particular timeevolution of the SALI allow us to distinguish fast and safely between regular and chaotic motion. Nevertheless, we have to define a specific numerical threshold value for determining the transition from order to chaos. After conducting extensive numerical experiments, integrating many sets of orbits, we conclude that a safe threshold value for the SALI is the value . Thus, in order to decide whether an orbit is regular or chaotic, one may follow the usual method according to which we check after a certain and predefined time interval of numerical integration, if the value of SALI has become less than the established threshold value. Therefore, if SALI the orbit is chaotic, while if SALI the orbit is regular thus making the distinction between regular and chaotic motion clear and beyond any doubt. For the computation of SALI we used the LPVI
code [CMD14, ], a fully operational routine which efficiently computes a suite of many chaos indicators for dynamical systems in any number of dimensions.
In our computations, we set time units as a maximum time of numerical integration. The vast majority of orbits (regular and chaotic) however, need considerable less time to find one of the two exits in the limiting curve and eventually escape from the system (obviously, the numerical integration is effectively ended when an orbit passes through one of the escape channels and escapes). Nevertheless, we decided to use such a vast integration time just to be sure that all orbits have enough time in order to escape. Remember, that there are the so called “sticky orbits" which behave as regular ones during long periods of time. Here we should clarify, that orbits which do not escape after a numerical integration of time units ( yr) are considered as nonescaping or trapped. In fact, orbits with escape periods equal to many Hubble times are completely irrelevant to our investigation since they lack physical meaning.
4 Numerical results
Figure  Type  

3a  1:1 loop  4.10  0.00  0  0.01  100  
3a  1:1 loop  1.85  0.00  0  0.01  100  
3b  1:3 box  0.00  1.25  0  0.10  100  
3c  chaotic  1.00  0.00  0  0.01  164  
3d  chaotic  2.40  0.00  0  0.01  129 
The main objective of our investigation is to determine which orbits escape and which remain trapped, distinguishing simultaneously between regular and chaotic trapped motion^{3}^{3}3Generally, any dynamical method requires a sufficient time interval of numerical integration in order to distinguish safely between ordered and chaotic motion. Therefore, if the escape rate of an orbit is very low or even worse if the orbit escapes directly from the system then, any chaos indicator (the SALI in our case) will fail to work properly due to insufficient integration time. Hence, it is pointless to speak of regular or chaotic escaping orbits.. Furthermore, two additional properties of the orbits will be examined: (i) the channels through which the stars escape and (ii) the timescale of the escapes (we shall also use the terms escape period or escape rates). In the present paper, we shall explore these dynamical quantities for various values of the energy, always within the interval .
Our numerical calculations indicate that apart from the escaping orbits there is always a considerable amount of nonescaping orbits. In general terms, the majority of nonescaping regions corresponds to initial conditions of regular orbits, where a third integral of motion is present, restricting their accessible phase space and therefore hinders their escape. However, there are also chaotic orbits which do not escape within the predefined interval of time units and remain trapped for vast periods until they eventually escape to infinity. At this point, it should be emphasized and clarified that these trapped chaotic orbits cannot be considered, by no means, neither as sticky orbits nor as super sticky orbits with sticky periods larger than time units. Sticky orbits are those who behave regularly for long time periods before their true chaotic nature is fully revealed. In our case on the other hand, this type of orbits exhibit chaoticity very quickly as it takes no more than about 100 time units for the SALI to cross the threshold value (SALI ), thus identifying beyond any doubt their chaotic character. Therefore, we decided to classify the initial conditions of orbits in both the physical and phase space into four main categories: (i) orbits that escape through , (ii) orbits that escape through , (iii) nonescaping regular orbits and (iv) trapped chaotic orbits.
Additional numerical computations reveal that the nonescaping regular orbits are mainly 1:1 loop orbits for which a third integral applies^{4}^{4}4The total angular momentum is an approximately conserved quantity (integral of motion), even for orbits in non spherical potentials., while other types of secondary resonant orbits are also present. In Fig. 3a we present two 1:1 thin loop orbits, while in Fig. 3b a typical example of a secondary 1:3 resonant quasiperiodic orbit is shown. The notation we use for the regular orbits is according to [CA98, , ZC13, ], where the ratio of those integers corresponds to the ratio of the main frequencies of the orbit, where main frequency is the frequency of greatest amplitude in each coordinate. Main amplitudes, when having a rational ratio, define the resonances of an orbit. Finally in Figs. 3cd we observe two orbits escaping through (exit channel 1) and (exit channel 2), respectively. All regular orbits shown in Figs. 3ab were computed until time units, while on the other hand, the escaping orbits presented in Figs. 3cd were calculated for 5 time units more than the corresponding escape period in order to visualize better the escape trail. In Table 1 we provide the type, the exact initial conditions and the value of the energy for all the depicted orbits.
A simple qualitative way for distinguishing between regular and chaotic motion in Hamiltonian systems is by plotting the successive intersections of the orbits using a Poincaré Surface of Section (PSS) [HH64, ]. In particular, a PSS is a twodimensional (2D) slice of the entire fourdimensional (4D) phase space. The following Fig. 4(ab) shows two PSSs at the critical Jacobi energy . Fig. 4a corresponds to the physical space showing orbits crossing with , while orbits crossing with at the phase space are presented in Fig. 4b. In both cases, the initial conditions of the orbits are integrated forwards in time plotting a dot at each crossing through the surface of section. The outermost black thick curve circumscribing each PSS is the limiting curve which in the physical space is defined as , while in the phase space is given by
(10) 
It is seen, that the vast majority of the surfaces of section is covered by a unified and extended chaotic sea however, stability regions with initial conditions corresponding to regular orbits are also present. In these stability islands the quasiperiodic orbits admit a third integral of motion, also known as “adelphic integral" (see e.g., [C63, ]), which hinders the test particles (stars) from escaping. Moreover, we observe that some areas in the chaotic sea of the plane are less densely occupied than others. We remark that in the following grid exploration in the physical and phase space, we do not consider all initial conditions in the same chaotic domain as one the same orbit.
4.1 The structure of physical space
Our exploration begins in the physical space and in the left column of Fig. 5 we present the orbital structure of the plane for three values of the energy, where the initial conditions of the orbits are classified into four categories by using different colors. Specifically, gray color corresponds to regular nonescaping orbits, blue color corresponds to trapped chaotic orbits, green color corresponds to orbits escaping through channel 1, while the initial conditions of orbits escaping through exit channel 2 are marked with red color. The vertical, dashed, magenta lines indicate the position of the two Lagrangian points and . For , that is an energy level just above the critical escape energy , we see that the vast majority of initial conditions corresponds to escaping orbits however, a stability island of nonescaping regular orbits is present at the central upper part of the interior galactic region. We also see that the entire interior region is completely fractal which means that there is a highly dependence of the escape mechanism on the particular initial conditions of the orbits. In other words, a minor change in the initial conditions has as a result the star to escape through the opposite exit channel, which is of course, a classical indication of chaotic motion. With a much closer look at the plane we can identify some additional tiny stability islands which are deeply buried in the escape domain. As we proceed to higher energy levels the fractal area reduces and several basins of escape start to emerge. By the term basin of escape, we refer to a local set of initial conditions that corresponds to a certain escape channel. Indeed for we observe the presence of escape basins which mainly have the shape of thin elongated bands. Moreover, the extent of the central stability island seems to be unaffected by the increase in the orbital energy. With increasing energy the interior region in the physical spaces becomes less and less fractal and broad, welldefined basins of escape dominate when . The fractal regions on the other hand, are confined mainly near the boundaries between the escape basins or in the vicinity of the stability islands. Furthermore, it is seen that for apart from the central stability island, there is also a set of three smaller islands of regular motion corresponding to 1:3 resonant orbits. It should be emphasized that our classification shows that in all three energy levels trapped chaotic motion is almost negligible as the initial conditions of such orbits appear only as lonely points around the boundaries of the stability islands. Here we must point out three interesting phenomena that take place with increasing energy: (i) the regions of forbidden motion are significantly reduced, (ii) the two escape channels become more and more wide and (iii) the amount of orbits with initial conditions outside the interior region (especially outside ) that escape from the opposite exit channel increases.
In the right column of Fig. 5 we show how the escape times of orbits are distributed on the physical space. The escape time is defined as the time when a star passes one of the two vertical, dashed lines at with velocity pointing outwards. Light reddish colors correspond to fast escaping orbits with short escape periods, dark blue/purpe colors indicate large escape rates, while white color denote both nonescaping regular and trapped chaotic orbits. The scale on the color bar is logarithmic. It is evident, that orbits with initial conditions close to the boundaries of the stability islands need significant amount of time in order to escape from the galaxy, while on the other hand, inside the basins of escape where there is no dependence on the initial conditions whatsoever, we measured the shortest escape rates of the orbits. We observe that for the escape periods of orbits with initial conditions in the fractal region are huge corresponding to tens of thousands of time units. As the value of the total orbital energy increases however, the escape times of orbits reduce significantly. In fact for and the basins of escape can be distinguished in the grids of Fig. 5, being the regions with reddish colors indicating extremely fast escaping orbits. Our numerical calculations indicate that orbits with initial conditions inside the basins have significantly small escape periods of less than 10 time units.
The following Fig. 6 presents the evolution of the percentages of the different types of orbits on the physical space when the dimensionless energy parameter varies. At this point we must clarify that we decided to merge the percentages of nonescaping regular and trapped chaotic orbits together because our computations indicate that always the rate of trapped chaotic orbits is extremely small (less than 0.1%) and therefore it does not contribute to the overall orbital structure of the dynamical system. One may observe, that when , that is just above the critical escape energy , escaping orbits through and share about 96% of the available physical space, so the two exit channels can be considered equiprobable. As the value of the energy increases however, the percentages of escaping orbits through exit channels 1 and 2 start to diverge, especially for . In particular, the amount of orbits that escape through exhibits a linear increase, while on the other hand, the rate of escaping orbits through the opposite exit channel displays a linear decrease. At the highest energy level studied escaping orbits through is the most populated family occupying about 55% of the physical plane, while only about 40% of the same plane corresponds to initial conditions of orbits that escape through . Furthermore, the percentage of nonescaping (regular plus chaotic) orbits is much less affected by the shifting of the energy displaying a rather constant value around 5%. Thus taking into account all the abovementioned analysis we may conclude that at low energy levels where the fractility of the physical space is maximum the stars do not show any particular preference regarding the escape channel, while on the contrary, at high enough energy levels where basins of escape dominate it seems that exit channel 2 is more preferable.
4.2 The structure of the phase space
We continue our investigation in the phase space and we follow the same numerical approach as discussed previously. In Fig. 7 we depict the orbital structure of the plane for three values of the energy, using different colors in order to distinguish between the four main types of orbits (nonescaping regular; trapped chaotic; escaping through and escaping through ). Just above the escape energy, that is when , we see that the structure of the phase plane is very similar to that of the corresponding physical plane shown in Fig. 5. It is observed that about 90% of the plane is occupied by initial conditions of escaping orbits. Moreover, the vast escape domain is highly fractal, while basins of escape, if any, are negligible. In all studied cases the areas of regular motion correspond mainly to retrograde orbits (i.e., when a star revolves around the galaxy in the opposite sense with respect to the motion of the galaxy itself), while there are also some smaller stability islands of prograde orbits at high energy levels (i.e., for ). The area on the phase plane covered by escape basins grows drastically with increasing energy and at high enough energy levels the fractal regions are confined mainly at the boundaries of the stability islands of nonescaping orbits. In particular, the escape basins first appear at mediocre values of energy inside the fractal, extended escape region having the shape of thin elongated bands however, as we proceed to higher energy levels they grow in size dominating in the phase space as wellformed broad domains. It is worth noticing the initial conditions of trapped chaotic orbits at the central region of the stability islands for which indicates the presence of a thin chaotic layer (separatrix^{5}^{5}5The separatrix prevents chaotic orbits from escaping in the same manner that an adelphic integral hinders quasiperiodic orbits from escaping.). The Coriolis forces dictated by the rotation of the compact active nucleus of the Seyfert galaxy makes the phase planes to be asymmetric with respect to the axis and this phenomenon is usually known as “Coriolis asymmetry" (e.g., [I80, ]).
The distribution of the corresponding escape times of orbits on the phase space as a function of the energy parameter is shown in the right column Fig. 7. One may observe that the results are very similar to those presented earlier in Fig. 5, where we found that orbits with initial conditions inside the basins of escape have the smallest escape rates, while on the other hand, the longest escape times correspond to orbits with initial conditions either in the fractal regions of the plots, or near the boundaries of the islands of regular motion. Another interesting way of measuring the escape period of an orbit is by counting how many intersection the orbit has with the axis before it escapes. The regions of the phase space for in Fig. 8 are colored according to the number of intersections with the axis upwards . Our computations reveal that orbits that escape directly without any intersection with the axis (0 intersections) or performing between 3 and 10 intersections or more than 10 intersections share about three fourths of the escape region. Orbits that cross one time the axis before escape correspond to about 15% of the escape realm, while only about 10% of the escaping orbits perform two intersections before leaving the system passing through one of the Lagrangian points. It is evident, that orbits with initial conditions located near the boundaries of the stability islands perform numerous intersections with the axis before they eventually escape to infinity. On the contrary, orbits with initial conditions in the escape basins need only a couple of intersection until they escape.
The evolution of the percentages of the different types of orbits on the phase space as a function of the dimensionless energy parameter is presented in Fig. 9. Once more, we have to point out that the percentages of nonescaping ordered and trapped chaotic orbits are merged together in a single trend line because the percentage of trapped chaotic orbits is extremely small (less than 0.5%) throughout. It is seen, that at the lowest energy level studied escaping orbits through and share about 90% of the phase plane, while only 10% of the same plane corresponds to initial conditions of nonescaping (regular plus chaotic) orbits. Furthermore, we observe that for the rates of all types of orbits remain practically unperturbed by the shifting on the value of the orbital energy. For larger values of the energy however, the percentages of escaping orbits exhibit a small divergence, while the portion of nonescaping orbits displays a minor increase and for they reach about 12.5%. Specifically, the rate of orbits escaping through slightly reduces, while that of orbits escaping through continues the monotone behavior. It should be emphasized that the divergence of the percentages of escaping orbits observed in the phase space is much smaller than that found to exist in the physical space. Therefore, one may reasonably deduce that in the phase space the increase in the value of the energy does not influence significantly the orbital content of the dynamical system, thus leaving the percentages of the orbits almost the same throughout. In addition, since that the rates of escaping orbits remain unperturbed within the energy range, we may say that the two exit channels can be considered almost equiprobable in the phase space.
4.3 An overview analysis
The colorcoded grids in the physical as well as in the phase space provide sufficient information on the phase space mixing however, for only a fixed value of the Jacobi constant. Hénon back in the late 60s [H69, ], introduced a new type of plane which can provide information not only about stability and chaotic regions but also about areas of trapped and escaping orbits using the section , (see also [BBS08, ]). In other words, all the orbits of the stars of the galaxy are launched from the axis with , parallel to the axis . Consequently, in contrast to the previously discussed types of planes, only orbits with pericenters on the axis are included and therefore, the value of the dimensionless energy parameter can be used as an ordinate. In this way, we can monitor how the energy influences the overall orbital structure of our dynamical system using a continuous spectrum of energy values rather than few discrete energy levels. In Fig. 10a we present the orbital structure of the plane when , while in Fig. 10c the distribution of the corresponding escape times of orbits is depicted.
We observe that for relatively low energy values the interior galactic region is highly fractal, while basins of escape are located only at the outer parts of the plane, that is outside and (or in other words in the exterior galactic region). Furthermore, we can identify the presence of two small stability islands of prograde nonescaping regular orbits which were invisible in the general plots shown earlier in Figs. 5 and 7. For larger values of energy however, the structure of the plane changes drastically and the most important differences are the following: (i) several basins of escape are formed inside the fractal escape region, (ii) apart from the main stability island a second smaller one emerges near . It should be pointed out that in the blowups of the diagram several additional extremely tiny islands of stability have been identified^{6}^{6}6An infinite number of regions of (stable) quasiperiodic (or small scale chaotic) motion is expected from classical chaos theory., (iii) at high enough energy levels we see that a thin chaotic layer composed of initial conditions of trapped chaotic orbits appear at the central region of the main stability island and (iv) the extent of the main island of regular motion grows rapidly and for the position of initial conditions of nonescaping regular orbits exceed the interior region . The following Table 2 shows the percentages of the four main types of orbits in the colorcoded grids shown in Fig. 10(ab).
Figure  Trapped  Nonescaping  Channel 1  Channel 2 

10a  0.34  19.05  31.94  48.67 
10b  1.46  28.72  26.53  43.29 
In order to obtain a more complete view on how the total orbital energy influences the nature of orbits in our galaxy model, we follow a similar numerical approach to that explained before but now all orbits of stars are initiated from the vertical axis with . In particular, this time we use the section , , launching orbits parallel to the axis. This allow us to construct again a 2D plane in which the coordinate of orbits is the abscissa, while the logarithmic value of the energy is the ordinate. The orbital structure of the plane when is shown in Fig. 10b. The vertical, dashed, magenta lines indicate the position of the Lagrangian points and . The black solid line is the limiting curve which distinguishes between regions of allowed and forbidden motion and is defined as
(11) 
A very complicated orbital structure is reveled in the plane which however, in general terms, is very similar with that of the plane. One may observe that for the forbidden regions of motion around and completely disappear. In fact the value is another critical energy level. It is seen that above that critical value the left part of the plane is occupied almost completely with initial conditions of orbits that escape through , while the right part of the same plane on the other hand, contains mostly initial conditions of nonescaping regular orbits.
It is evident from the results presented in Figs. 10(cd) that the escape times of the orbits are strongly correlated to the escape basins. In addition, one may conclude that the smallest escape periods correspond to orbits with initial conditions inside the escape basins, while orbits initiated in the fractal regions of the planes or near the boundaries of stability islands have the highest escape rates. In both types of planes the escape times of orbits are significantly reduced with increasing energy. Combining all the numerical outcomes presented in Figs. 10(ad) we may say that the key factor that determines and controls the escape times of the orbits is the value of the orbital energy (the higher the energy level the shorter the escape rates), while the fractality of the basin boundaries varies strongly both as a function of the energy and of the spatial variable.
The evolution of the average value of the escape time of orbits as a function of the dimensionless energy parameter is given in Fig. 11a. It is seen, that for low values of energy the average escape time of orbits reduces rapidly, while for , it seems to saturate around 80 time units. We feel it is important to justify this behaviour of the escape time. As previously explained, a Lyapunov orbit acts like a barrier, reaches the ZVC on both sides of the opening and returns along the same path thus, connecting two opposite branches of the ZVC. Therefore, we could exploit this fact in order to quantify the escape process. For this purpose, we define as the distance between the two points at which the Lyapunov orbit touches the two opposite branches of the ZVC. Using we can measure, in a way, the width of the escape channels. Fig. 11b depicts the evolution of as a function of the energy parameter . We observe an almost exponential growth of with increasing energy. This means that as the value of the Jacobi constant increases the escape channels (which are of course symmetrical) become more and more wide and therefore, orbits need less and less time in order to find one of the two openings in the open equipotential surface and eventually escape from the system. This geometrical feature explains why for low values of energy orbits consume large time periods wandering inside the open equipotential surface until they eventually locate one of the two exits and escape to infinity. Therefore, we may conclude that it is the geometrical shape of the exit channels that justifies the decreasing pattern of shown in Fig. 11a.
5 Formation of spiral arms
It is well known that when stars escape from star clusters through the Lagrangian points form complicated structures known as tidal tails or tidal arms (e.g., [CMM05, , dMCM05, , JBPE09, , KMH08, , KKBH10, ]). In the same vein, in the case of barred spiral galaxies we observe the formation of spiral arms (e.g., [EP14, , GKC12, , MT97, , MQ06, , QDB11, , RDQ08, , SK91, ]). Usually a star cluster rotates around its parent galaxy in a circular orbit thus, the tidal tails follow the curvature defined by the circular orbit and they are nearly straights for orbits moving in large galactocentric distances (see e.g., [JBPE09, ]). The spiral arms of barred galaxies on the other hand, are formed by the nonaxisymmetric perturbation of the bar. In particular, in barred spiral galaxies with prominent spirallike shape the two arms start from the ends of the bar and then they wind up around the bananashaped forbidden regions which enclose the Lagrangian points and .
Almost all Seyfert galaxies are spiral galaxies and have been among the most intensively studied objects in astronomy (e.g., [SSN14, ]). Therefore, it would be very challenging to investigate if our simple gravitational model has the ability to realistically simulate the formation of spiral arms. As a first step, we integrated the two orbits shown in Fig. 3(cd) for 7 time units after crossing the Lagrangian points and in Fig. 12 we present the corresponding paths of the orbits that escape from opposite exits. It is seen, that both orbits follow indeed a spiral path around the yellow bananashaped isolines of the effective potential in which the Lagrangian points and are enclosed.
Taking into account that our initial tests are positive regarding the formation of the spiral structure, we reconfigured our integration routine for the computation of the basins of escape to yield output of all orbit trajectories for a threedimensional grid of size for . This dense grid of initial conditions is centered at the origin of the phase space and we allow for all orbits both sings of . Our purpose is to monitor the timeevolution of the path of the orbits and simulate the spiral arm formation. At all orbits are regularly distributed within the Lagrangian radius inside the interior galactic region. In Fig. 13 we present four time snapshots of the position of the stars in the physical plane. We observe that at time units, which corresponds to 125 Myr the majority of stars are still inside the interior galactic region however, a small portion of stars has already escaped passing through either or . At time units (250 Myr) we have the first indication of the formation of spiral arms, where with green and red color we depict stars that escaped through exit channels 1 and 2, respectively, while stars that remain inside the Lagrangian radius are shown in gray. As time goes by we see that the two symmetrical spiral arms grow in size and at time units (500 Myr) the Seyfert galaxy has obtained its complete spiral structure, while the interior galactic region is uniformly filled. At this point we have to point out that the standard integration code, which was used to obtain the results presented in Figs. 5 and 7, plots a point at every step of the numerical integration. In Fig. 13(ad) on the other hand, the density of points along one star orbit is taken to be proportional to the velocity of the star. Being more specific, a point is plotted (showing the position of a star), if an integer counter variable which is increased by one at every integration step, exceeds the velocity of the star. Following this technique we can simulate, in a way, a real body simulation of the spiral evolution of the galaxy, where the density of stars will be highest where the corresponding velocity is lowest. Therefore, we proved that for the scenario of evolution of spiral arms in our galaxy model is indeed viable. Additional numerical simulations reveal that this is also true for other (lower or higher) values of energy. In fact, the particular energy level mainly controls how tight the spiral arms wind up around the forbidden regions of motion.
In section 2 we explained that the values of the mass and the scale length entering Eq. (1) were chosen having in mind the Seyfert galaxy NGC 1566 with total mass equal to 1.2 [KKB05, ]. Here we have to point out that NGC 1566 is not just any Seyfert galaxy; it is the second brightest Seyfert galaxy known. It is also the brightest and most dominant member of the Dorado Group, a loose concentration of galaxies located approximately 40 million lightyears away that together comprise one of the richest galaxy groups of the southern hemisphere. NGC 1566 is an intermediate spiral galaxy of type SAB(rs)bc, meaning that while it does not have a well defined barshaped region of stars at its central region  like normal barred spirals  it is not quite an unbarred spiral either. The small but extremely bright nucleus of NGC 1566 is clearly visible in Fig. 14 [img, ], a telltale sign of its membership of the Seyfert class of galaxies. The centres of such galaxies are very active and luminous, emitting strong bursts of radiation and potentially harbouring supermassive black holes that are many millions of times the mass of the Sun. This image highlights the beauty and aweinspiring nature of this unique galaxy group, with NGC 1566 glittering and glowing, its bright nucleus framed by swirling and symmetrical lavender arms. In Fig. 14 we have superposed the time snapshot at time units (500 Myr) of our simulation above the real image of the galaxy. We observe that the central interior region adequately matches the active galactic core, while the two spiral arms obtained from the numerical simulation fit the corresponding real ones. Thus, we may conclude that our simple onecomponent gravitational effective potential (2) can, in a way, model not only the escape dynamics but also the spiral formation and evolution of Seyfert galaxies.
Undoubtedly, a Seyfert galaxy is a very complex stellar entity and therefore, we need to assume some necessary simplifications and assumptions in order to be able to study mathematically the orbital behavior of such a complicated stellar system. For this purpose, our gravitational model is intentionally simple and contrived, in order to give us the ability to study all the different aspects regarding the kinematics and dynamics of the system. Nevertheless, contrived models can surely provide an insight into more realistic stellar systems, which unfortunately are very difficult to be explored, if we take into account all the astrophysical aspects (i.e., gas, mergers, close encounters etc). Selfconsistent models on the other hand, are mainly used when conducting body simulations. However, this application is entirely out of the scope of the present paper. Once more, we have to point out that the simplicity of our model is necessary, otherwise it would be extremely difficult, or even impossible, to apply the extensive and detailed numerical calculations presented in this study. Our singlecomponent potential neglect many features which must exist in real Seyfert galaxies such as colors/ages of stars populations, deviations from simple symmetries and timeindependent perturbations. However, it is still significant that the conclusions derived from this simple potential all agree, both qualitatively and semiquantitatively. This suggests that our model is robust, depending only on such topological properties of the phase space.
6 Discussion and conclusions
The aim of this work was to shed some light to the trapped or escaping nature of orbits in a dynamical gravitational model describing a Seyfert galaxy. In order to describe the properties of the Seyfert galaxy as a whole we used an asymmetric Plummertype potential. At this point we should emphasize that a Seyfert galaxy is composed of a disk, bulge a central nucleus and a bar. However, our dynamical system has two degrees of freedom and for this type of Hamiltonians we usually adopt single component potentials to model the properties of barred galaxies (e.g., [BT08, , C00, , CK98, , C81, , C83, , CP80, ]) or galaxies in general (e.g., [BBP06, , CP03a, , CP03b, , CM85, , KC01, , KC02, , MAB95, , VWD12, ]). This approach allow us to directly determine the effects of the single component potential to the the different types (resonances) of orbits, the amount of chaos, the escape properties and other dynamical aspects.
The Jacobi integral of motion has a critical threshold value (or escape energy) which determines if the escape process is possible or not. In particular, for energies smaller than the threshold value, the equipotential surface is closed and it is certainly true that escape is impossible. For energy levels larger than the escape energy however, the equipotential surface opens and two exit channels appear through which the particles can escape to infinity. Here it should be emphasized, that if a test particle does has energy larger than the escape value, there is no guarantee that the star will certainly escape from the system and even if escape does occur, the time required for an orbit to transit through an exit and hence escape to infinity may be vary long compared with the natural escaping time. We managed to distinguish between ordered/chaotic and trapped/escaping orbits and we also located the basins of escape leading to different exit channels, finding correlations with the corresponding escape times of the orbits. Our extensive and thorough numerical investigation strongly suggests, that the overall escape mechanism is a very complicated procedure and very dependent on the value of the Jacobi integral.
Since a distribution function of the model was not available so as to use it for extracting the different samples of orbits, we had to follow an alternative path. We defined for several values of the Jacobi integral, dense grids of initial conditions and regularly distributed in the area allowed by the energy level on the physical and phase space, respectively and then we identified regions of order/chaos and bound/escape. In each type of grid the step separation of the initial conditions along the axes (or in other words the density of the grid) was controlled in such a way that there are always at least orbits to be integrated and classified. For the numerical integration of the orbits in each type of grid, we needed about between 1 hour and 5 days of CPU time on a Pentium DualCore 2.2 GHz PC, depending on the escape rates of orbits in each case. For each initial condition, the maximum time of the numerical integration was set to be equal to time units however, when a particle escaped the numerical integration was effectively ended and proceeded to the next available initial condition.
The present article provides quantitative information regarding the escape dynamics in spiral Seyfert galaxies. The main numerical results of our research can be summarized as follows:

In all examined cases, areas of bounded motion and regions of initial conditions leading to escape in a given direction (basins of escape), were found to exist in both the physical and the phase space. The several escape basins are very intricately interwoven and they appear either as welldefined broad regions or thin elongated bands. Regions of bounded orbits first and foremost correspond to stability islands of regular orbits where a third integral of motion is present.

A strong correlation between the extent of the basins of escape and the value of the Jacobi integral was found to exist. Indeed, for low energy levels the structure of both the physical and the phase space exhibits a large degree of fractalization and therefore the majority of orbits escape choosing randomly exit channels. As the value of the energy increases however, the structure becomes less and less fractal and several basins of escape emerge. The extent of these basins of escape is more prominent at high energy levels.

In both the physical and the phase space we identified a small portion of trapped chaotic orbits which do not escape within the predefined maximum time of numerical integration. The initial conditions of these orbits are located either in thin chaotic layers (separatrix) inside the regular areas of motion, or near the boundaries of stability islands. Moreover, these orbits reveal their chaotic nature very quickly, so they cannot be considered as sticky orbits.

We observed, that in many cases the escape process is highly sensitive dependent on the initial conditions, which means that a minor change in the initial conditions of an orbit lead the test particle to escape through the other exit channel. These regions are the exact opposite of the escape basins, are completely intertwined with respect to each other (fractal structure) and are mainly located in the vicinity of stability islands. This sensitivity towards slight changes in the initial conditions in the fractal regions implies, that it is impossible to predict through which exit the particle will escape.

Our calculations revealed, that the escape times of orbits are directly linked to the basins of escape. In particular, inside the basins of escape as well as relatively away from the fractal domains, the shortest escape rates of the orbits had been measured. On the other hand, the longest escape periods correspond to initial conditions of orbits either near the boundaries between the escape basins or in the vicinity of the stability islands.

It was detected, that for energy levels slightly above the escape energy the majority of the escaping orbits have considerable long escape rates (or escape periods), while as we proceed to higher energies the proportion of fast escaping orbits increases significantly. This behavior was justified, by taken into account a geometrical feature of the equipotential surface. Specifically, for low values of energy the width of the exit channels is too small therefore, stellar orbits consume large time periods wandering inside the open equipotential surface until they eventually locate one of the two exits and escape to infinity. Moreover, it was proved that the width of the escape channels grows rapidly with increasing energy, while at the same time the average escape time decreases significantly.

We presented evidence that escaping stars with values of energy higher than the critical Jacobi value form spiral arms outside the interior galactic region defined by the Lagrangian radius. This procedure has a striking similarity to the formation of tidal tails in star clusters. We also conducted numerical simulations proving that our gravitational effective potential can model the formation as well as the evolution of spiral arms in real Seyfert galaxies.
Judging by the detailed and novel outcomes we may say that our task has been successfully completed. We hope that the present numerical analysis and the corresponding results to be useful in the field of escape dynamics in Active Galactic Nuclei (AGN). The outcomes as well as the conclusions of the present research are considered, as an initial effort and also as a promising step in the task of understanding the escape mechanism of orbits in AGNs. Taking into account that our results are encouraging, it is in our future plans to modify properly our galactic model in order to expand our investigation into three dimensions and explore the entire sixdimensional phase space. In addition, body simulations may elaborate the formation as well as the evolution of spiral structure in Seyfert galaxies.
Acknowledgments
The author would like to thank the two anonymous referees for the careful reading of the manuscript and for all the apt suggestions and comments which allowed us to improve both the quality and the clarity of the paper.
References
 (1) Aguirre, J., Vallego, J.C., Sanjuán, M.A.F.: Wada basins and chaotic invariant sets in the HénonHeiles system. Phys. Rev. E 64, 0662081–11 (2001)
 (2) Aguirre, J., Sanjuán, M.A.F.: Limit of small exits in open Hamiltonian systems. Phys. Rev. E 67, 0562011–7 (2003)
 (3) Aguirre, J., Vallego, J.C., Sanjuán, M.A.F.: Wada basins and unpredictability in Hamiltonian and dissipative systems. Int. J. Mod. Phys. B 17, 41714175 (2003)
 (4) Aguirre, J., Viana, R.L., Sanjuán, M.A.F.: Fractal structures in nonlinear dynamics. Rev. Mod. Phys. 81, 333386 (2009)
 (5) Aguirregabiria, J.M.: Chaotic scattering around black holes. Phys. Lett. A 224, 234238 (1997)
 (6) Altmann, E.G., Portela, J.S.E., Tél, T.: Leaking chaotic systems. Rev. Mod. Phys. 85, 869918 (2013)
 (7) Barrio, R., Blesa, F., Serrano, S.: Fractal structures in the HénonHeiles Hamiltonian. Europhys. Lett. 82, 100031–6 (2008)
 (8) Barrio, R., Blesa, F., Serrano, S.: Bifurcations and safe regions in open Hamiltonians. New J. Phys. 11, 0530041–12 (2009)
 (9) Belmonte, C., Boccaletti, D., Pucacco, G.: Stability of axial orbits in galactic potentials. Celest. Mech. Dyn. Astron. 95, 101116 (2006)
 (10) Benet, L., Trautman, D., Seligman, T.: Chaotic scattering in the Restricted ThreeBody Problem. I. The Copenhagen Problem. Celest. Mech. Dyn. Astron. 66, 203228 (1996)
 (11) Benet, L., Seligman, T., Trautman, D.: Chaotic scattering in the Restricted ThreeBody Problem II. Small mass parameters. Celest. Mech. Dyn. Astron. 71, 167189 (1998)
 (12) Binney, J., Tremaine, S.: Galactic Dynamics, Princeton Univ. Press, Princeton, USA (2008)
 (13) Bleher, S., Grebogi, C., Ott, E., Brown, R.: Fractal boundaries for exit in Hamiltonian dynamics. Phys. Rev. A 38, 930938 (1988)
 (14) Bleher, S., Grebogi, C., Ott, E.: Bifurcation to chaotic scattering. Phys. D 46, 87121 (1990)
 (15) Bleher, S., Ott, E., Grebogi, C.: Routes to chaotic scattering. Phys. Rev. Let. 63, 919922 (1989)
 (16) Boyd, P.T., McMillan, S.L.W.: Initialvalue space structure in irregular gravitational scattering. Phys. Rev. A 46, 62776287 (1992)
 (17) Capuzzo Dolcetta, R., Di Matteo, P., Miocchi, P.: Formation and evolution of clumpy tidal tails around globular clusters. Astron. J. 129, 19061921 (2005)
 (18) Caranicolas, N.D.: Maps describing motion in strong bars. New Astronomy 5, 397402 (2000)
 (19) Caranicolas, N.D., Karanis, G.I.: Chaos in barred galaxy models. Astrophys. Space Sci. 259, 4556 (1998)
 (20) Caranicolas, N.D., Papadopoulos, N.J.: Connecting gravitational potential parameters to chaos in elliptical galaxies. New Astronomy 9, 103110 (2003)
 (21) Caranicolas, N.D., Papadopoulos, N.J.: Comparing maps to symplectic integrators in a galactictype Hamiltonian. J. Astrophys. Astr. 24, 8597 (2003)
 (22) Caranicolas, N.D., Papadopoulos, N.J.: The S(c) spectrum machine to visualize the motion in galaxies. Astron. Nachr. 328, 556561 (2007)
 (23) Carpintero, D.D., Aguilar, L.A.: Orbit classification in arbitrary 2D and 3D potentials. Mon. Not. R. Astron. Soc. 298, 121 (1998)
 (24) Carpintero, D.D., Maffione, N., Darriba, L.: LPVI code: A program to compute a suite of variational chaos indicators. Astronomy and Computing 5, 1927 (2014)
 (25) Chen, Q., Ding, M., Ott, E.: Chaotic scattering in several dimensions. Phys. Lett. A 145, 93100 (1990)
 (26) Churchill, R.C., et al. in Como Conference Proceedings on Stochastic Behavior in Classical and Quantum Hamiltonian Systems, Volume 93, Lecture Notes in Physics, ed. G. Casati, J. Fords (Berlin: Springer), 76 (1979)
 (27) Contopoulos, G.: On the existence of a third integral of motion. Astron. J. 68, 114 (1963)
 (28) Contopoulos, G.: Asymptotic curves and escapes in Hamiltonian systems. Astron. Astrophys. 231, 4155 (1990)
 (29) Contopoulos, G.: The effects of resonances near coroation in barred galaxies. Astron. Astrophys. 102, 265278 (1981)
 (30) Contopoulos, G.: Bifurcations, gaps and stochasticity in barred galaxies. Astrophys. J. 275, 511528 (1983)
 (31) Contopoulos, G.: Order and Chaos in Dynamical Astronomy. Springer, Berlin (2002)
 (32) Contopoulos, G., Kaufmann, D.: Types of escapes in a simple Hamiltonian system. Astron. Astrophys. 253, 379388 (1992)
 (33) Contopoulos, G., Kandrup, H.E., Kaufmann, D.: Fractal properties of escape from a twodimensional potential. Phys. D 64, 310323 (1993)
 (34) Contopoulos, G., Magnenat, P.: Simple threedimensional periodic orbits in a galactictype potential. Cel. Mech. 37, 387414 (1985)
 (35) Contopoulos, G., Harsoula, M., LukesGerakopoulos, G.: Periodic orbits and escapes in dynamical systems. Celest. Mech. Dyn. Astron. 113, 255278 (2012)
 (36) Contopoulos, G., Papayannopoulos, Th.: Orbits in weak and strong bars. Astron. Astrophys. 92, 3346 (1980)
 (37) de Moura, A.P.S., Letelier, P.S.: Fractal basins in HénonHeiles and other polynomial potentials. Phys. Lett. A 256, 362368 (1999)
 (38) Di Matteo, P., Capuzzo Dolcetta, R., Miocchi, P.: Clumpy substructures in globular cluster tidal tails. Celest. Mech. Dyn. Astron. 91, 5973 (2005)
 (39) Ding, M., Grebogi, C., Ott, E., Yorke, J.A.: Transition to chaotic scattering. Phys. Rev. A 42, 70257040 (1990)
 (40) Eckhardt, B.: Fractal properties of scattering singularities. J. Phys. A 20, 59715979 (1987)
 (41) Eckhardt, B.: Irregular scattering. Phys. D 33, 8998 (1988)
 (42) Eckhardt, B., Jung, C.: Regular and irregular potential scattering. J. Phys. A 19, L829L833 (1986)
 (43) Ernst, A., Just, A., Spurzem, R., Porth, O.: Escape from the vicinity of fractal basin boundaries of a star cluster. Mon. Not. R. Astron. Soc. 383, 897906 (2008)
 (44) Ernst, A., Peters, T.: Fractal basins of escape and the formation of spiral arms in a galactic potential with a bar. Mon. Not. R. Astron. Soc. 443, 25792589 (2014)
 (45) Gaspard, P., Rice, S.A.: Scattering from a classically chaotic repellor. J. Chem. Phys. 90, 22252241 (1989)
 (46) Grand, R.J.J., Kawata, D., Cropper, M.: The dynamics of stars around spiral arms. Mon. Not. R. Astron. Soc. 421, 15291538 (2012)
 (47) Hénon, M.: Numerical exploration of the restricted problem, V. Astron. Astrophys. 1, 223238 (1969)
 (48) Hénon, M.: Chaotic scattering modelled by an inclined billiard. Phys. D 33, 132156 (1988)
 (49) Hénon, M., Heiles, C.: The applicability of the third integral of motion: Some numerical experiments. Astron. J. 69, 7379 (1964)
 (50) Hut, P.: The topology of threebody scattering. Astron. J. 88, 15491559 (1983)
 (51) Hut, P., Bahcall, J.N.: Binarysingle star scattering. I  Numerical experiments for equal masses. Astrophys. J. 268, 319341 (1983)
 (52) Innanen K.A.: The Coriolis asymmetry in the classical restricted 3body problem and the Jacobian integral. Astron. J. 85, 8185 (1980)
 (53) José, J.V., Rojas, C., Saletan, E.J.: Elastic particle scattering from two hard disks. Amer. J. Phys. 60, 587592 (1992)
 (54) Jung, C.: Can the integrability of Hamiltonian systems be decided by the knowledge of scattering data? J. Phys. A 20, 17191732 (1987)
 (55) Jung, C., Lipp, C., Seligman, T.H.: The inverse scattering problem for chaotic Hamiltonian systems. Ann. Phys. 275, 151189 (1999)
 (56) Jung, C., MejiaMonasterio, C., Seligman, T.H.: Scattering one step from chaos. Phys. Lett. A 198, 306314 (1995)
 (57) Jung, C., Pott, S.: Classical cross section for chaotic potential scattering. J. Phys. A 22, 29252938 (1989)
 (58) Jung, C., Richter, P.H.: Classical chaotic scatteringperiodic orbits, symmetries, multifractal invariant sets. J. Phys. A 23, 28472866 (1990)
 (59) Jung, C., Scholz, H.J.: Cantor set structures in the singularities of classical potential scattering. J. Phys. A 20, 36073618 (1987)
 (60) Jung, C., Tel, T.: Dimension and escape rate of chaotic scattering from classical and semiclassical cross section data. J. Phys. A 24, 27932805 (1991)
 (61) Just, A., Berczik, P., Petrov, M., Ernst, A.: Quantitative analysis of clumps in the tidal tails of star clusters. Mon. Not. R. Astron. Soc. 392, 969981 (2009)
 (62) Kandrup, H.E., Siopis, C., Contopoulos, G., Dvorak, R.: Diffusion and scaling in escapes from twodegreesoffreedom Hamiltonian systems. Chaos 9, 381392 (1999)
 (63) Karanis, G.I., Caranicolas, N.D.: Transition from regular motion to chaos in a logarithmic potential. Astron. Astrophys. 367, 443448 (2001)
 (64) Karanis, G.I., Caranicolas, N.D.: A new dynamical spectrum for galactic potentials. Astron. Nachr. 323, 311 (2002)
 (65) Kennedy, J., Yorke, J.A.: Basins of Wada. Phys. D 51, 213225 (1991)
 (66) Kilborn, V.A., Koribalski, B.S., Forbes, D.A., Barnes, D.G., Musgrave, R.C.: A widefield HI study of the NGC 1566 group. Mon. Not. R. Astron. Soc. 356, 7788 (2005)
 (67) Küpper, A.H.W., Macleod, A., Heggie, D.C.: On the structure of tidal tails. Mon. Not. R. Astron. Soc. 387, 12481252 (2008)
 (68) Küpper, A.H.W., Kroupa, P., Baumgardt, H., Heggie, D.C.: Tidal tails of star clusters. Mon. Not. R. Astron. Soc. 401, 105120 (2010)
 (69) Lai, Y.C., de Moura, A.P.S., Grebogi, C.: Topology of highdimensional chaotic scattering. Phys. Rev. E 62, 64216428 (2000)
 (70) Lai, Y.C., Grebogi, C., Blümel, R., Kan, I.: Crisis in chaotic scattering. Phys. Rev. Let. 71, 22122215 (1993)
 (71) Lau, Y.T., Finn, J.M., Ott, E.: Fractal dimension in nonhyperbolic chaotic scattering. Phys. Rev. Let. 66, 978981 (1991)
 (72) Lipp, C., Jung, C.: From scattering singularities to the partition of a horseshoe. Chaos 9, 706714 (1999)
 (73) Mahon, M.E., Abernathy, R.A., Bradley, B.O., Kandrup, H.E.: Transient ensemble dynamics in timeindependent galactic potentials. Mon. Not. R. Astron. Soc. 275, 443453 (1995)
 (74) Masset, F., Tagger, M.: Nonlinear coupling of spiral waves in disk galaxies: a numerical study. Astron. Astrophys. 322, 442454 (1997)
 (75) Minchev, I., Quillen, A.C.: Radial heating of a galactic disc by multiple spiral density waves. Mon. Not. R. Astron. Soc. 368, 623636 (2006)
 (76) Moser, J.: On the generalisation of a theorem of A. Liapunoff. Commun. Pure Appl. Math. 11, 257271 (1958)
 (77) Motter, A.E., Lai, Y.C.: Dissipative chaotic scattering. Phys. Rev. E 65, 015205 (2002)
 (78) Navarro, J.F., Henrard, J.: Spiral windows for escaping stars. Astron. Astrophys. 369, 11121121 (2001)
 (79) Nagler, J.: Crash test for the Copenhagen problem. Phys. Rev. E 69, 066218 (2004)
 (80) Nagler, J.: Crash test for the restricted threebody problem. Phys. Rev. E 71, 026227 (2005)
 (81) Ott, E.: Chaos in dynamical systems, 2nd Ed., Cambridge University Press (2002)
 (82) Papadopoulos, N.J., Caranicolas, N.D.: Do active galaxies have a massive halo component? New Astron. 12, 1119 (2006)
 (83) Poon, L., Campos, J., Ott, E., Grebogi, C.: Wada basins boundaries in chaotic scattering. Int. J. Bifurc. Chaos 6, 251266 (1996)
 (84) Press, H.P., Teukolsky, S.A, Vetterling, W.T., Flannery, B.P.: Numerical Recipes in FORTRAN 77, 2nd Ed., Cambridge Univ. Press, Cambridge, USA (1992)
 (85) Quillen, A.C., Dougherty, J., Bagley, M.B., Minchev, I., Comparetta, J.: Structure in phase space associated with spiral and bar density waves in an Nbody hybrid galactic disc. Mon. Not. R. Astron. Soc. 417, 762784 (2011)
 (86) Roškar, R., Debattista, V.P., Quinn, T.R., Stinson, G.S., Wadsley, J.: Riding the spiral waves: Implications of stellar migration for the properties of galactic disks. Astrophys. J. 684, L79L82 (2008)
 (87) Rückerl, B., Jung, C.: Scaling properties of a scattering system with an incomplete horseshoe. J. Phys. A 27, 5577 (1994)
 (88) Schneider, J., Tél, T., Neufeld, Z.: Dynamics of “leaking” Hamiltonian systems. Phys. Rev. E 66, 0662181–6 (2002)
 (89) SchnorrMüller, A., StorchiBergmann, T., Nagar, N.M., Ferrari, F.: Gas inflows towards the nucleus of the active galaxy NGC 7213. Mon. Not. R. Astron. Soc. 438, 33223331 (2014)
 (90) Sellwood, J.A., Kahn, F.D.: Spiral modes driven by narrow features in angularmomentum density. Mon. Not. R. Astron. Soc. 250, 278299 (1991)
 (91) Seoane, J.M., Aguirre, J., Sanjuán, M.A.F., Lai, Y.C.: Basin topology in dissipative chaotic scattering. Chaos 16, 0231011–8 (2006)
 (92) Seoane, J.M., Sanjuán, M.A.F., Lai, Y.C.: Fractal dimension in dissipative chaotic scattering. Phys. Rev. E 76, 0162081–6 (2007)
 (93) Seoane, J.M., Sanjuán, M.A.F.: Exponential decay and scaling laws in noisy chaotic scattering. Phys. Lett. A 372, 110116 (2008)
 (94) Seoane, J.M., Huang, L., Sanjuán, M.A.F., Lai, Y.C.: Effects of noise on chaotic scattering. Phys. Rev. E 79, 0472021–4 (2009)
 (95) Seoane, J.M., Sanjuán, M.A.F.: Escaping dynamics in the presence of dissipation and noisy in scattering systems. Int. J. Bifurc. Chaos 9, 27832793 (2010)
 (96) Seoane, J.M., Sanjuán, M.A.F.: New developments in classical chaotic scattering. Rep. Prog. Phys. 76, 016001 (2013)
 (97) Simó C.: Dynamical properties in Hamiltonian Systems. Applications to Celestial Mechanics. Text of the lectures delivered at the Centre de Recerca Matemàtica on January 2731, (2014)
 (98) Siopis, C.V., Contopoulos, G., Kandrup, H.E.: Escape probabilities in a Hamiltonian with two channels of escape. New York Acad. Sci. Ann. 751, 205212 (1995)
 (99) Siopis, C.V., Kandrup, H.E., Contopoulos, G., Dvorak, R.: Universal properties of escape. New York Acad. Sci. Ann. 773, 221230 (1995)
 (100) Siopis, C.V., Kandrup, H.E., Contopoulos, G., Dvorak, R.: Universal properties of escape in dynamical systems. Celest. Mech. Dyn. Astron. 65, 57681 (1996)
 (101) Skokos, C.: Alignment indices: A new, simple method for determining the ordered or chaotic nature of orbits. J. Phys. A: Math. Gen. 34, 1002910043 (2001)
 (102) Sweet, D., Ott, E.: Fractal basin boundaries in higherdimensional chaotic scattering. Phys. Lett. A 266, 134139 (2000)
 (103) Szebehely, V.: Theory of Orbits. Academic Press, New York (1967)
 (104) Terra, M.O., de Assis, S.C.: Escape dynamics and fractal basin boundaries in the planar Earth–Moon system. Celest. Mech. Dyn. Astron. 120, 105130 (2014)
 (105) The Hubble European Space Agency Information Centre: http://www.spacetelescope.org/
 (106) Valluri, S.R., Wiegert, P.A., Drozd, J., Da Silva, M.: A study of the orbits of the logarithmic potential for galaxies. Mon. Not. R. Astron. Soc. 427, 23922400 (2012)
 (107) Zotos, E.E.: Trapped and escaping orbits in axially symmetric galactictype potential. PASA 29, 161173 (2012)
 (108) Zotos, E.E.: Order and chaos in a galactic model with a strong nuclear bar. Research in Astr. Astrophys. 12, 500512 (2012)
 (109) Zotos, E.E.: Revealing the evolution, the stability and the escapes of families of resonant periodic orbits in Hamiltonian systems. Nonlin. Dyn. 73, 931962 (2013)
 (110) Zotos, E.E.: A Hamiltonian system of three degrees of freedom with eight channels of escape: The Great Escape. Nonlin. Dyn. 76, 13011326 (2014)
 (111) Zotos, E.E.: Escapes in Hamiltonian systems with multiple exit channels: Part I Nonlin. Dyn. 78, 13891420 (2014)
 (112) Zotos, E.E., Carpintero, D.D.: Orbit classification in the meridional plane of a disk galaxy model with a spherical nucleus. Celest. Mech. Dyn. Astron. 116, 417438 (2013)