© The Authors 2025

1. Introduction

In optical images, bars are observed in almost half of all nearby disk galaxies (see, e.g. Reese et al. 2007; Barazza et al. 2008). This percentage increases to nearly 70% when near-infrared images are considered (Knapen et al. 2000; Menéndez-Delmestre et al. 2007).

Structures in barred galaxies have to be supported by stellar orbits (Contopoulos 2002; Binney & Tremaine 2008). The orbital study of barred galaxies has been the subject of numerous papers. It is well known that 2D and 3D regular orbits connected to the x₁ stable periodic family of orbits (Contopoulos & Papayannopoulos 1980) support the shape of the bar in face-on as well as edge-on projections (e.g. Patsis et al. 2003, 2009; Patsis & Katsanikas 2014; Chaves-Velasquez et al. 2017) in barred galaxies.

On the other hand, many studies in recent years have shown that the envelope of the bar, as well as the ring and spiral structures outside the bar, can be supported by sticky chaotic orbits that behave like regular ones for many dynamical times of the system (Kaufmann & Contopoulos 1996; Patsis et al. 1997; Tsoutsis et al. 2008, 2009; Harsoula & Kalapotharakos 2009; Contopoulos & Harsoula 2013; Tsigaridi & Patsis 2015). These chaotic orbits are sticky along unstable manifolds of unstable periodic orbits near the region of corotation in the case of 2D galactic models. These studies were carried out within the framework of the manifold theory introduced by Voglis et al. (2006) and Romero-Gómez et al. (2006), which is the result of the outflow of sticky chaotic orbits connected with the unstable dynamics around the bar’s unstable Lagrangian points, L₁ (L₂) (see also Danby 1965).

In models of barred galaxies, corotation is usually placed close to the end of the bar (Contopoulos 1980; Athanassoula 1980). Bars are classified according to their ratio, R, between the corotation radius and the bar length, into ultrafast bars (R < 1.0), fast bars (1.0 < R < 1.4), and slow bars (R > 1.4).

In this paper we study the role of 2D and 3D sticky chaotic orbits¹ in a Milky Way-like barred galactic model of three degrees of freedom, outside the region of corotation, in supporting features like ring structures for a relatively small perturbation of the bar. These orbits are related to the planar Lyapunov orbits, PL₁ (PL₂) (which are the continuation of the unstable Lagrangian equilibrium points, L_1, 2, for energy levels above E_j, L1, 2), or to other 3D regular or simply unstable periodic orbits, all located outside corotation for Jacobi integrals, E_j, above E_j, L1, 2.

We study the case of a relatively weak bar with a maximum perturbation of ≈12% in forces (see Section 3) and a constant pattern speed, Ω_bar = 45 km/s/kpc, which lies in the range of values proposed for the Milky Way by Gerhard (2011) and Sanders et al. (2019). We show that 2D sticky chaotic orbits related to the planar PL₁ (PL₂) family, as well as other unstable periodic orbits in the region of corotation, support ring and spiral structures, while 3D small perturbations of these orbits support only ring structures of R₁ or R₁R₂ types (Buta & Combes 1996) for long enough times compared to the Hubble time. These rings can be considered as ‘outer rings’ and are placed near the outer Lindblad resonance (OLR) (Athanassoula & Bosma 1985; Buta & Combes 1996). The outer rings are closely linked to the spiral pattern and are occasionally wrapped pairs of spiral arms that nearly close together. In such cases, the ring is referred to as an outer pseudoring and is represented by R′ (Buta & Combes 1996).

We use the 4D Poincaré surface of section (PSS) of the 6D phase space to investigate the different types of regular or sticky chaotic orbits of the barred galactic model. Directly visualizing the phase space of Hamiltonian systems with three degrees of freedom is challenging because their PSS is four dimensional. Since Froeschlé’s pioneering work (Froeschlé 1970), various methods have been developed to facilitate the visualization of 4D PSSs. Common approaches include 2D and 3D projections of the 4D PSS (Martinet & Magnenat 1981; Vrahatis et al. 1996, 1997), stereoscopic representations (Martinet et al. 1981; Richter et al. 2014), and Lagrangian descriptors (Agaoglou et al. 2021). In this study we use the method of colour and rotation (Patsis & Zachilas 1994). This method was used to study various structures of the phase space of galactic-type Hamiltonian systems with three degrees of freedom (Katsanikas & Patsis 2011, 2022; Katsanikas et al. 2011a,b, 2013; Moges et al. 2024).

Particularly important measures for the study of sticky chaotic orbits are the time evolution of their finite-time Lyapunov characteristic number (FT-LCN) and the type and speed of diffusion compared to the dynamical time of the barred galaxy or the Hubble time. Such a comparison shows that this outwards diffusion is very slow and that these sticky chaotic orbits behave in practice like regular ones for many dynamical times of the system supporting the morphological structures of the galaxy. After this stickiness time, these orbits get quickly diffused outside the system through superdiffusing procedures. This underlines the importance of stickiness phenomena along unstable asymptotic manifolds or other hypersurfaces in the 6D phase space outside corotation.

In Sect. 2 we describe the galactic model of a Milky Way-like barred galaxy. In Sect. 3 we give a quantitative estimate of the force strength of the bar. In Section 4 we make an exploration of the 6D phase space of the model with the help of 4D PSSs. We study the 3D sticky chaotic orbits that support ring-like structures by plotting their apocentric sections on the plane of rotation. We also study the time evolution of the FT − LCN and the type of diffusion of these sticky chaotic orbits. We present our conclusions in Sect. 5.

2. The model

2.1. The axisymmetric potential

The mass distribution of the stellar disk of a barred galaxy is very well represented by an exponential profile (Freeman 1970). The profile of the surface density, Σ(R), of an exponential disk has the form

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{0}exp(-R/R_d),\end{array}$$(1)

where Σ₀ is the central surface density and R_d is the disk’s scale length. The vertical density distribution used in stellar disks usually has a sech form and therefore the space density distribution of the galactic disk has the form

$$\begin{array}{c}\hfill \rho (R,z)={\rho }_{0}exp(-R/R_d){\mathrm{sech}}^n(-|z|/z_0),\end{array}$$(2)

where z₀ is the scale height and n is typically ≈1–3.

The potential that corresponds to an exponential density profile of form (2) is the well-known Miyamoto-Nagai (MN) disk (Miyamoto & Nagai 1975)

$$\begin{array}{c}\hfill V_d(R,z)=\frac{-G{M}_d}{\sqrt{R^2+{(a_d+\sqrt{z^2+b_d^2})}^2}},\end{array}$$(3)

where M_d is the total mass of the disk, a_d is the radial scale length, and b_d is the vertical scale length (Binney & Tremaine 2008, p. 73–74). For a_d = 0 the potential (3) becomes a spherical Plummer distribution, while for b_d = 0 it becomes an infinitely thin Kuzmin disk (Kuzmin 1956). The parameters a_d and b_d control the range of thickness of the mass distribution.

However, the surface density corresponding to a single MN potential is too low near the centre and too high at large radii (e.g. Smith et al. 2015). A more realistic axisymmetric potential for the stellar disk of a barred galaxy is a combination of two MN disks for the thin and the thick parts of the galactic disk (e.g. Katsanikas & Patsis 2011):

$$\begin{array}{cc}\hfill V_{d_1(R,z)}& =\frac{-G{M}_{d_1}}{\sqrt{R^2+{(a_{d_1}+\sqrt{z^2+b_{d_1}^2})}^2}}\hfill \\ \hfill V_{d_2(R,z)}& =\frac{-G{M}_{d_2}}{\sqrt{R^2+{(a_{d_2}+\sqrt{z^2+b_{d_2}^2})}^2}},\hfill \end{array}$$(4)

where R = x² + y² is the radius on the plane of rotation. The potential

$$\begin{array}{c}\hfill V_{\mathrm{ax}}=V_{d_1(R,z)}+V_{d_2(R,z)}\end{array}$$(5)

can be considered as a good approximation of the axisymmetric part of the Milky Way potential. The parameters in Eq. (4) have the following values: M_d1 = 2.05 × 10¹⁰ M_⊙, M_d2 = 25.47 × 10¹⁰ M_⊙a_d1 = 0 kpc, a_d2 = 7.258 kpc, b_d1 = 0.495 kpc, and b_d2 = 0.520 kpc. These values were proposed by Miyamoto & Nagai (1975) in order to fit the corresponding rotation curve to the Galactic one.

2.2. The bar potential

The bar potential is a Milky Way-like bar potential proposed by Long & Murali (1992):

$$\begin{array}{c}\hfill V_b=\frac{G{M}_b}{2a}ln\left(\frac{y-a+T_{-}}{y+a+T_{+}}\right),\end{array}$$(6)

with $T_{\pm }=\sqrt{[{(a\pm y)}^2+x^2+{(b+\sqrt{c^2+z^2})}^2]}$, a = 5.25 kpc (the length of the major axis of the bar on the plane of rotation, along the y-axis), b = 2.1 kpc (the length of the minor axis of the bar on the plane of rotation, along the x-axis), and c = 1.6 kpc (the thickness of the bar along the z-axis). The length of the bar is within the range of values given in the literature for the Milky Way’s bar (Wegg et al. 2015; Lucey et al. 2023). The ratios between the lengths of the bar’s axes are consistent with the ones given in the literature (Gerhard 2002; Rattenbury et al. 2007; Cao et al. 2013). M_b is the mass of the bar in 10¹⁰ M_⊙.

This potential has been used in the literature in N-body simulations and models of Milky Way-like galaxies (Lee et al. 1999; Pettitt et al. 2014; Efthymiopoulos et al. 2020). The length unit R = 1 corresponds to 1 kpc. The unit of mass M = 1 corresponds to 10¹⁰ M_⊙. By setting the Newton’s constant equal to unity, G = 1, the time unit t_un = 1 corresponds to 1.49 × 10¹⁴ sec ≈ 5M year. The Hubble time is ≈3000t_un. Finally, the unit of velocity is calculated as V_un = 207.4 km/sec.

3. Force strength of the bar

In order to give a measure of the bar’s perturbation we used the Q_f-force strength commonly used in the literature. The definition of the Q_f-force strength is given (Laurikainen et al. 2004) as

$$\begin{array}{c}\hfill Q_f(r)=\frac{|F_T{(r,\varphi )|}_{max}}{|F_R(r)|},\end{array}$$(7)

where $|F_R(r)|=\frac{\partial V_{\mathrm{ax}}}{\partial r}$ denotes the axisymmetric force at each radius on the x − y plane and |F_T(r, ϕ)|_max is the maximum value azimuthally of the tangential force $|F_T(r,\varphi )|=(1/r)\frac{\partial V_{\mathrm{ax}}}{\partial \varphi }$ due to the perturbation of the bar at each radius of the x − y plane.

In Figure 1 we plot the value of Q_f as a function of the radius, r, for M_b = 5 × 10¹⁰ M_⊙ (the mass of the bar) in Eq. (6). The maximum value of Q_f(r) corresponds to a perturbation of 12%, which lies within the range of the corresponding values observed in real barred galaxies. The range of the Q-force perturbation as derived from observations of real barred galaxies varies approximately between the values of 10% and 70% (see Laurikainen et al. 2004 and Buta et al. 2005, 2009). Consequently, the value of 12% that we adopt here as our Q-force perturbation is a rather small value for a bar’s strength. The vertical line corresponds to the position of the L₁ Lagrangian equilibrium point that lies shortly after the end of the bar.

Fig. 1. Q-force perturbation of the bar, Qf, as a function of the radius, r, of the bar. The maximum Qf is approximately 12%. The vertical line gives the position of the Lagrangian points L1, 2 that are located shortly after the end of the bar.

4. Exploration of the phase space

4.1. Locating the Lagrangian points L₁, L₂

We fixed the pattern speed of the bar to Ω_b = 45 km/s/kpc in real units and Ω_b = 0.217 in the units of the system described above. The Hamiltonian of the system, which rotates around the z-axis with angular velocity Ω_b, can be written in Cartesian coordinates in the form

$$\begin{array}{c}\hfill H=\frac{p_x^2}{2}+\frac{p_y^2}{2}+\frac{p_z^2}{2}+V(x,y,z)-{\mathrm{\Omega }}_b(x{p}_y-y{p}_x),\end{array}$$(8)

where V(x, y, z) = V_ax + V_b and p_x, p_y, and p_z are the canonically conjugate momenta. Hereafter, we will denote the numerical value of the Hamiltonian by E_j and refer to it as the Jacobi constant or, more loosely, as the ‘energy’. The corresponding equations of motion (Hamilton equations) are

$$\begin{array}{cc}\hfill \dot{x}& =\frac{\mathit{dx}}{\mathit{dt}}=p_x+{\mathrm{\Omega }}_{b}y\hfill \\ \hfill \dot{y}& =\frac{\mathit{dy}}{\mathit{dt}}=p_y-{\mathrm{\Omega }}_{b}x\hfill \\ \hfill \dot{z}& =\frac{\mathit{dz}}{\mathit{dt}}=p_z\hfill \\ \hfill \dot{p_x}& =\frac{d{p}_x}{\mathit{dt}}=-\frac{\vartheta V}{\vartheta x}+{\mathrm{\Omega }}_b{p}_y\hfill \\ \hfill \dot{p_y}& =\frac{d{p}_y}{\mathit{dt}}=-\frac{\vartheta V}{\vartheta y}-{\mathrm{\Omega }}_b{p}_x\hfill \\ \hfill \dot{p_z}& =\frac{d{p}_z}{\mathit{dt}}=-\frac{\vartheta V}{\vartheta z}·\hfill \end{array}$$(9)

The Poincaré surface of section for a specific Jacobi constant E_j is four dimensional (4D). The equations of motion are solved for a given E_j, with initial conditions $(x_0=0,\dot{x_0}<0,y_0,{\dot{y}}_0,z_0,{\dot{z}}_0)$ in the rotating frame of reference (with the help of Eqs. (9)).

We first located the Lagrangian equilibrium points, L_1, 2, as well as the Jacobi constant, E_L1, 2, that corresponds to these points. For this purpose we used the equations of motion (9) and solved the following system of equations:

$$\begin{array}{c}\hfill \frac{\mathit{dx}}{\mathit{dt}}=\frac{\mathit{dy}}{\mathit{dt}}=\frac{\mathit{dz}}{\mathit{dt}}=\frac{d{p}_x}{\mathit{dt}}=\frac{d{p}_y}{\mathit{dt}}=\frac{d{p}_z}{\mathit{dt}}=0.\end{array}$$(10)

Equation (10) gives five different solutions that correspond to the five Lagrangian equilibrium points: two unstable (L₁,  L₂) located close to the end of the bar along its major axis, one stable (L₃) at the centre of the galaxy, and two stable (L₄,  L₅) located at a perpendicular direction to the main axis of the bar. By selecting convenient initial conditions as guess values we find the solution for the L₁ Lagrangian equilibrium point: x_L1 = 0,  y_L1 ≈ 6.3,  z_L1 = 0,  p_xL1 ≈ −1.37,  p_yL1 = 0,  p_zL1 = 0. The corresponding Jacobi constant was calculated from Eq. (8): E_jL1 = −4.55622. The Lagrangian equilibrium point L₂ is symmetric to the L₁ point with respect to the x-axis.

If we consider the radius of the Lagrangian point L₁ close to the corotation radius, then the ratio of the corotation radius to the length of the bar (taken as equal to its major axis, a = 5.25, as defined in Section 2.2) is R ≈ 1.2 < 1.4 and this can be considered as a fast bar. In a study of 225 real observed barred galaxies by Geron et al. (2023) a classification of bars is made according to their ratio, R, between the corotation radius and the bar length, into ultrafast bars (R < 1.0; 11% of their sample), fast bars (1.0 < R < 1.4; 27%), and slow bars (R > 1.4; 62%). In contrast, (Sierra et al. 2015) find results consistent with the values of R obtained from the computer simulations of Rautiainen et al. (2008), indicating that 34 out of 57 galaxies (60%) have fast bars.

4.2. 4D Poincaré surfaces of section

We now explore the 6D phase space for Jacobi constants, E_j, close to and above E_jL1 = −4.55622. We chose the values E_j = −4.5, −4.0, −3.5.

According to the Lyapunov subcentre theorem (Rabinowitz 1982; Weinstein 1973; Moser 1976) there exist at least two families of periodic orbits that are the continuation of the Lagrangian equilibrium points L_1, 2 for energies above E_jL1, 2. These families are i) the planar family of periodic orbits named PL_1, 2 and ii) the vertical family of periodic orbits named vertical Lyapunov orbits, VLO_1, 2. For this work we studied the planar periodic orbits, PL_1, 2, for energies above E_jL1, 2. A study of the vertical Lyapunov orbits, VLO_1, 2, has been undertaken recently (Katsanikas & Patsis 2022).

We study the linear stability of the planar periodic orbits, PL_1, 2, with the help of the variational equations considering small deviations from the initial conditions $(y_{P{L}_1},\dot{y_{P{L}_1}},x_{P{L}_1}=0,\dot{x_{P{L}_1}}<0)$ (see Appendix A). The PL_1, 2 family consists of two branches, denoted branch 1 and branch 2. The corresponding stability indices for these branches are (Fig. 2a) (a) the vertical stability indices, b1(1) (branch 1) and b1(2) (branch 2) (black curves) corresponding to the stability of planar periodic orbits with respect to perturbations perpendicular to the orbital plane and (b) the planar stability indices, b2(1) (branch 1) and b2(2) (branch 2) (red curves) corresponding to the stability with respect to perturbations on the orbital plane. We plotted the stability indices as a function of the energy, E_n (Fig. 2a). These branches arise due to a saddle-node bifurcation, as indicated by an arrow in Fig. 2a. This bifurcation is a higher-dimensional extension of the saddle-node bifurcation commonly found in Hamiltonian systems with two degrees of freedom.

Fig. 2. (a) Stability indices, b1, b2, of the planar PL1, 2 orbit (branch 1 and branch 2) as a function of the Jacobi constant (or energy), En, for bar mass Mb = 5. The vertical stability index, b1, (black) corresponds to the stability of planar periodic orbits with respect to perturbations perpendicular to the orbital plane, while the planar stability index, b2, (red) corresponds to the stability with respect to perturbations on the orbital plane. The three vertical lines correspond to the energy levels at which we investigated the phase space. (b) Characteristics of the two planar branches of the PL1 family that are joined and disappear at the bifurcation point. (c) Branch 1 and (d) branch 2 planar PL1 orbits plotted with the same colour as the corresponding energies (vertical lines) in Fig. 2a.

We give the conditions of the stability of the planar periodic orbits here (for more details see Appendix A). The periodic orbit is stable if b₁ < 0, b₂ < 0, simply unstable if b₁ < 0, b₂ > 0 or b₁ > 0, b₂ < 0, double unstable if b₁ > 0, b₂ > 0 and complex unstable if b₁ and b₂ are complex numbers.

In Hamiltonian systems with two degrees of freedom, a saddle-node bifurcation occurs when the following two conditions coexist: Firstly, in the stability diagram, which represents the stability index of a periodic orbit as a function of the energy (or another non-linearity parameter), a family of periodic orbits that is initially unstable (branch 1) coexists with a family of periodic orbits that is initially stable (branch 2) and in a specific value of energy, the two branches merge and disappear (the bifurcation point). Secondly, the bifurcation point coincides with the maximum or minimum of the characteristic curve of the two branches. The characteristic curve describes the position of periodic orbits in a particular direction (such as the y-direction in our case) within the PSS as a function of the energy.

In our case, where the Hamiltonian is of three degrees of freedom, the same conditions must be satisfied for the planar PL_1, 2 family of orbits. The first condition described above holds here, but instead of stability and instability (as seen in systems with two degrees of freedom), we observe initially that branch 1 of the PL_1, 2 family is simply unstable (also referred to as semi-instability by Broucke 1969), while branch 2 is stable. Finally, the two branches merge and disappear at the bifurcation point marked by an arrow in Fig. 2a. The branch 1 family approaches the bifurcation point as simply unstable, while the branch 2 family does so as double unstable. The second condition described above is also satisfied, as the bifurcation point coincides with the minimum of the characteristic curve of branch 1 and the maximum of the characteristic curve of branch 2 (Fig 2b).

Following the stability of branch 1, in Fig. 2a, from smaller to larger values of energy, we observe that it is initially simply unstable, as expected, with −2 < b1(1) < 2 and b2(1) <  − 2, as this family is the continuation of the unstable equilibrium points L_1, 2 for larger values of energy. Then it becomes stable at the point where the b2(1) index crosses the axis b2 = −2 and then simply unstable again when the stability index b1(1) crosses the axis b1 = −2, just before the bifurcation point indicated by an arrow.

Branch 2 is initially stable and then transitions to simple instability when b1(2) crosses the axis b1 = 2 upwards and to stability again when b1(2) crosses the axis b1 = 2 downwards. Then, for greater energy values, branch 2 becomes, successively, simply unstable when b2(2) crosses the axis b2 = −2 downwards and double unstable when b1(2) crosses the axis b1 = −2 downwards. In Figure 2c, the planar PL₁ orbits of branch 1, and in Fig. 2d, the planar PL₁ orbits of branch 2, are plotted with the same colour as the corresponding energies (vertical lines) in Fig. 2a.

If we fix the initial values $(y_0,\dot{y_0},z_0,\dot{z_0})$ on the 4D PSS x = 0, the value of the coordinate ẋ₀ (velocity in the rotating frame) is calculated for a specific energy, E_n = E_j, with the help of Eqs. (8), (9) as

$$\begin{array}{c}\hfill \dot{x_0}=\pm \sqrt{2(E_n-V(x_0,y_0,z_0))-{\dot{y_0}}^2-{\dot{z_0}}^2+{\mathrm{\Omega }}_b^2(x_0^2+y_0^2)}.\end{array}$$(11)

We chose the negative solution of ẋ₀ in Eq. (11) for negative y₀, which is compatible with the direction of rotation of the orbits outside corotation, in the rotating frame of reference.

In Figure 3a,b we plot the 3D projection ($y,\dot{y},z$) of the 4D PSS and use the colour gradient to visualize the fourth dimension, ż, (using the method of colour and rotation introduced by Patsis & Zachilas 1994) of the 6D phase space. We plot in black the 2D PSS (y,ẏ, for x = 0 and $\dot{x}<0$) of the planar orbits that stay continuously on the x − y plane of rotation on the configuration space. The corresponding energy is E_j = −4.5.

Fig. 3. 4D PSS ($y,\dot{y},z,\dot{z}$) for x = 0 and $\dot{x}<0$, for energy level Ej = −4.5. (a) For smaller radii. (b) For larger radii. The colour gradient allows the visualization of the fourth dimension, ż, in rotational tori of regular orbits according to the colour bar on the left (minimum ż corresponds to 0 and maximum ż to 1). The rotational torus of the regular ring-like orbit is also shown. The 2D PSS of the planar orbits is plotted in black and the unstable manifolds of the planar unstable periodic orbits in the chaotic region are in red. 3D sticky chaotic orbits (green) escape from the system after very long times. (a′) and (b′) are the 2D projections (y,ẏ) of the 4D PSS of the same orbits as in (a) and (b).

The colour gradient enables the fourth dimension, ż, to be visualized in rotational tori of regular orbits (in Figs. 3a,b) according to the colour bar on the left (minimum ż corresponds to 0 and maximum ż to 1). We plot the 2D PSS of the planar orbits in black, and in red the unstable manifolds of the planar unstable periodic orbits in the chaotic region, which can escape to infinity. 3D sticky chaotic orbits (green in Fig. 3a) can escape from the system, but after a very long time compared to the Hubble time. Figs. 3a′,b′ are the 2D projections (y,ẏ) of the 4D PSS of the same orbits as in Figs. 3a,b, for a better understanding. Fig. 4 is the same as Fig. 3 for energy levels E_j = −4.0 in (a) and (a′) and E_j = −3.5 in (b) and (b′).

Fig. 4. Same as in Fig. 3, but for energies (a), (a′) Ej = −4.0 and (b), (b′) Ej = −3.5.

In this paper, we examine sticky chaotic orbits that persist within a specific region of the phase space for a significant time interval. These orbits must be distinguished from partially chaotic orbits, which also remain confined to a particular region, but differ in nature. Partially chaotic orbits obey an additional isolating integral beyond energy and are constrained by regular orbits (see Muzzio 2017, 2018; Carpintero & Muzzio 2020). The latter category, however, falls outside the scope of this study and will not be discussed further.

In Figs. 3, 4 we plot the two branches of the planar PL₂ periodic orbits (which is convenient, as we study the PSS for negative values of the y coordinate), which are always symmetric to the PL₁ family and have the same stability indices and characteristic curves (Figs. 2a,b). Moreover, the shape of the orbits is the same as for the PL₁ orbits shown in Figs. 2c,d, but symmetric with respect to the y = 0 axis.

In Fig. 3a we show the two branches of the PL₂ family (named PL₂ for branch 1 and PL′₂ for branch 2). We observe that the PL₂ orbit is unstable, while the PL′₂ orbit is stable (see Fig. 2a). The unstable manifold of the planar PL₂ family is found outside corotation, but well inside the region of organized motion, and cannot escape from the system. The 4D coloured rotational torus corresponds to a regular orbit that supports a ring around the bar of the galaxy and we call it the ‘ring orbit’ (see below). The green points correspond to a 3D sticky chaotic orbit with initial conditions above the last Kolmogorov-Arnold-Moser (KAM) curve (before the chaotic region) of the 2D PSS of the planar orbits. Although this orbit has a small z component (z = 0.02) and theoretically can escape through this dimension, it stays confined to a small region of the phase space due to stickiness phenomena and escapes after a very long time of integration (compared to the dynamical time of the system), supporting ring-like structures (see Section 4.4). This is also obvious in the 2D projection of the 4D PSS (Fig. 3a′).

We observe that for large radii far beyond corotation (Figs. 3b,b′), we find a plethora of regular orbits that correspond to rotational tori on this figure. These tori cannot however be populated by real particles because, apart from the 3D ring-like orbit (−8 < y < −7), all other regular orbits are found far enough from the end of the bar.

Figures 4a,a′ correspond to E_j = −4.0. Here again, the planar PL₂ orbit is unstable, while the PL′₂ one is stable. Sticky 3D chaotic orbits above the planar unstable manifold of the PL₂ orbit (green points), remain located there for a very long time before escaping through the third dimension. They present a very slow diffusion outwards (see Section 4.3) and again support ring-like structures on the configuration space (see Section 4.4).

Finally, Figs. 4b,b′ correspond to E_j = −3.5. At this energy level, both the PL₂ and PL′₂ orbits have become stable (see Fig. 2a). We plot in red the planar unstable manifold of an unstable periodic orbit in the chaotic region that is able to escape from the system and in green the 3D sticky chaotic orbits having a small initial z component above the last KAM curve of the regular motion of the 2D PSS.

In Figure 5, top, we plot the 3D ring orbit (red) on the configuration space (x, y, z) that supports an R₁ ring-like structure around the bar of the galaxy and whose rotational torus on the 4D PSS is shown in Fig. 3. This is a 3D regular orbit, while the corresponding 2D ring orbit on the x − y plane (black) is an unstable periodic orbit. This is shown in Fig. 5, bottom, where we plot the PSS (y,ẏ, for x = 0 and $\dot{x}<0$) of the planar orbits. We plot the 1D unstable manifold of the PL₂ orbit and some other planar orbits in black, together with the projection of the rotational torus of the ring orbit (red). We observe that the projection of the regular ring orbit is placed above the corresponding planar unstable periodic orbit of the same family. This phenomenon is possible when the vertical stability index, b₁, of the planar ring orbit indicates stability and the planar stability index, b₂, indicates instability (see Appendix A).

Fig. 5. Top: 3D regular orbit that supports an R1 ring-like structure of the galaxy, labelled ‘Ring Orbit’ (red), together with the planar unstable periodic orbit (black). Bottom: 2D PSS (y,ẏ) for x = 0 and $\dot{x}<0$. The projection of the 3D ring-like orbit is plotted in red. The energy level is Ej = −4.5.

4.3. Diffusion of chaotic orbits

We wanted to define a measure for the diffusion of the 3D chaotic orbits outwards. For this purpose we needed to first define a dynamical time of the galactic potential. A commonly used measure for the dynamical time of a barred galaxy is the half mass crossing time, T_hm, given by the relation

$$\begin{array}{c}\hfill T_{\mathrm{hm}}=\sqrt{\frac{2R_{\mathrm{hm}}^3}{\mathit{GMg}}},\end{array}$$(12)

where R_hm is the half mass radius of the system and M_g is the total mass of the barred galaxy, which is the sum of the bar’s mass, M_b, and the masses of the two exponential MN disks, M_d1 and M_d2 (see Section 2). In order to find R_hm, we used the following procedure. From the Poisson equation ▿²V(x, y, z) =  − 4πGρ(x, y, z) we determined the volume density, ρ(x, y, z), that corresponds to the total potential, V(x, y, z), of the barred galaxy. Then we found the surface density, σ(0, y), along the major axis of the bar by integrating the volume density, ρ(x, y, z), along the z direction:

$$\begin{array}{c}\hfill \sigma (0,y)={\displaystyle {\int }_{-z}^z}\rho (y,z)dz,\end{array}$$(13)

where −10 kpc < z < 10 kpc are the limits of integration. From Eq. (13) we could find the mass, M(r) = σ(r)π[(r + dr/2)² − (r − dr/2)²], of every spherical ring between the radii r − dr/2 and r + dr/2 and then find the cumulative mass, M_c(r), as a function of the radius. From this function we finally determined the half mass radius, R_hm, of the system. We find that R_hm is approximately 8.67 kpc and therefore the half mass crossing time, T_hm, is approximately 6.3 time units. The Hubble time corresponds to approximately 3000 time units (see Section 2.2). Hence, the Hubble time is about 470 times the T_hm of the galaxy, which is a realistic number as the lifetime of a galaxy corresponds to a few hundred crossing times, in general.

In Figure 6, the 2D PSS (y,ẏ) for x = 0 and $\dot{x}<0$ of planar orbits is plotted, in black, for an energy level E_j = −4.0 (see Figs. 4a,a′), which is above the energy, E_jL1, of the Lagrangian points L₁, L₂ (see Section 4.1). The islands of stability plotted in different colours correspond to the planar regular orbits shown in Fig. 7, corresponding to the same colours as in Fig. 6. They all support planar ring-like structures on the configuration space.

Fig. 6. 2D PSS (y,ẏ) for x = 0 and $\dot{x}<0$ of planar orbits (with $z,\dot{z}=0$) for energy Ej = −4.0. The coloured islands of stability correspond to orbits shown in Fig. 7.

Fig. 7. Planar orbits corresponding to the islands of stability of the same colours as in Fig. 6. These all support ring-like structures.

Next we studied the diffusion for three different cases of 3D sticky chaotic orbits for the same energy level as in Fig. 6. In Figure 8, the (y,ẏ) projection of the 4D PSS ($y,\dot{y},z,\dot{z}$) for x = 0 and $\dot{x}<0$ for three 3D sticky chaotic orbits is plotted, with different colours, and superimposed on the 2D PSS (y,ẏ) in black. We plot the projection of the 4D PSS of a 3D orbit initially lying slightly above a 2D island of stability in red, while in green we plot the projection of the 4D PSS of a 3D chaotic orbit lying above the unstable manifold of the planar PL₂ orbit. Lastly, the blue orbit is the projection of the 4D PSS of a 3D chaotic orbit lying slightly above an unstable manifold of an orbit just outside the last KAM curve of the 2D PSS of the planar orbits (black). All three orbits have initial conditions with a small z component (z = 0.02).

Fig. 8. 2D PSS (y,ẏ, for x = 0 and $\dot{x}<0$) of planar orbits for energy level Ej = −4.0 (black). The (y,ẏ) projections of the 4D PSS of the first 20 000 iterations of three 3D sticky chaotic orbits having a small third dimension z = 0.02 above the plane of rotation are plotted in colour. The red one behaves like a regular orbit and the green one like a weakly chaotic orbit for very long times. The blue one is diffused rather quickly outwards.

The time evolution of these three orbits on the 4D PSS ($y,\dot{y},z,\dot{z}$; for x = 0 and $\dot{x}<0$) is shown in Fig. 9, top. In fact, we plot the 3D projection ($y,\dot{y},z$) of the 4D PSS. On the top left we plot the PSS of the corresponding orbits with the same colours as in Fig. 8 for the first 20 000 iterations. The red orbit shows a regular orbital behaviour and the green one a sticky chaotic behaviour and both are located in a limited region of the PSS. More specifically, the red orbit lies upon the cyan islands of stability, in Fig. 6, and have the form of a rotational torus, while the green orbit is a sticky chaotic orbit that lies upon the chaotic region corresponding to the PL₂ unstable periodic orbit (see Fig. 4a,a′). Finally the blue orbit has already been diffused outwards (escape) during this time of integration. In the middle panel of Fig. 9, top, we integrate the orbits for a much longer time and plot them for iterations between 20 000 and 50 000. We observe that the red and green orbits have now been diffused in the z direction but that they still have not escaped from the galaxy. In the top right panel we plot the same orbits for a total of 150 000 iterations and see that they have been diffused inwards and fill a certain region closer to the centre of the galaxy. They still have not escaped outwards. This number of iterations corresponds to hundreds of times the Hubble time. In order to have a better understanding of the time evolution of these sticky chaotic orbits described above, we plot in Fig. 9, bottom, the 2D projections of the 4D PSS for the same periods of time as in Fig. 9, top. We plot the 2D PSS of the planar orbits in black, while in colour we plot the projections of the sticky 3D chaotic orbits presented above. In Fig. 9, bottom, on the right, green and red points coexist and occupy the same area of the PSS (the same volume in the 4D phase space), although only red points are shown in the figure due to projection phenomena.

Fig. 9. Top: 3D projection of the 4D PSS ($y,\dot{y},z$) of the 3D chaotic orbits plotted in colour in Fig. 8. Top left: Orbits for the first 20 000 iterations. The red and green orbits behave as ordered ones for very long times, while the blue one is diffused out of the system more rapidly. Top middle: Orbits for iterations between 20 000 and 50 000. Top right: Orbits for the next 100 000 iterations. The time evolutions of the red and green sticky chaotic orbits present a very slow diffusion inwards. Bottom: Corresponding 2D projections of the 4D PSS for the same time periods. The 2D PSS of the planar orbits is shown in black, with the projections of the sticky 3D chaotic orbits presented above shown in colour.

By calculating the time evolution of the finite-time Lyapunov characteristic number, FT-LCN, of these sticky chaotic orbits,

$$\begin{array}{c}\hfill FT-LCN=\frac{1}{t}ln\frac{w(t)}{w(0)},\end{array}$$(14)

where w(0) is an initial deviation vector of the orbit and w(t) is the deviation vector after a time t (see Appendix B for details), we can determine the time scales during which these sticky chaotic orbits behave like regular ones. The limit of FT-LCN for t → ∞ is the well-known Lyapunov characteristic number (LCN).

In Figure 10, top, we plot the 4D PSS for energy level E_j = −4.0, as in Fig. 9, but here in orange or yellow (the gradient of the colour indicates the ż component), we plot a sticky chaotic orbit that acts as an approximate barrier between the orbits that get diffused inwards (green) and the ones that get diffused outwards (blue). In bottom, we plot the time evolution of the FT-LCN of the two chaotic orbits described above (blue and green) on a logarithmic scale, using the same colours as the orbits in Figs. 8, 9. The Hubble time corresponds to 3000 time units. We observe that the blue chaotic orbit is sticky for some time and then (before a Hubble time) gets diffused quickly outwards, manifesting its chaotic behaviour and converges to a positive value of LCN. The orbit then escapes to very large distances from the system and therefore we cease the calculation of the corresponding FT-LCN. In contrast, the sticky chaotic orbit (green) behaves like a regular one for very long time periods (thousands of Hubble times) before manifesting its chaotic nature with the convergence of its FT-LCN at a positive value when it becomes diffused inwards and spread in a larger volume closer to the centre of the galaxy (see Fig. 9, on the right). It is remarkable that the corresponding FT-LCN of the green orbit presents an abrupt rise when it begins to spread in a much larger volume closer to the centre of the system (green chaotic orbit of Fig. 9). This sticky orbit manifests its chaotic nature at about 3000 Hubble times. We expect that after an extremely long time (that is not possible to follow in our numerical calculations) this orbit will finally escape from the system and get diffused outwards, like the blue orbit. The red orbit in Fig. 9 has the same behaviour and the same time evolution of its FT-LCN as the green one.

Fig. 10. Top: 4D PSS for energy level Ej = −4.0. The same as in Fig. 9, but here in orange or yellow (the gradient of the colour indicates the ż component) showing a sticky chaotic orbit that acts as an approximate barrier between the orbits that get diffused inwards (green) and those that get diffused outwards (blue). Bottom: Time evolution of FT-LCN for the two sticky chaotic orbits of Fig. 9. A chaotic orbit (blue) diffused faster outwards and reaching very large distances, and a sticky chaotic orbit (green) that behaves like a regular one for very long times and finally converges to a positive LCN.

We conclude that all these sticky chaotic orbits in the vicinity of corotation behave like regular ones in the lifetime of the galaxy and support morphological structures of the barred galaxy outside corotation (see Section 4.4) for extremely long times compared to the Hubble time. It seems that physical barriers (in the 4D PSS) exist between the sticky chaotic orbits that get diffused inwards and the ones that get diffused outwards and escape from the system (orange and yellow orbit in Fig. 9, on the left).

In Figure 11 we plot the 3D projection of the 4D PSS of the three sticky chaotic orbits described above, together with the surface of zero velocity (SZV; derived from the equation En = V(x, y, z)−(1/2)Ω_sp²(x² + y²), for En = −4.0), which limits the allowed area of motion of the orbits in this energy level (seen from two different angles). From this figure we understand that the red chaotic orbit is limited inside the SZV and probably needs an extremely long time to escape from it through the open bottle neck of the SZV. The limit in the z direction inside the SZV defines the thickness of the thick disk, which is consistent with the numbers given for our own Milky Way (Gilmore & Reid 1983). On the other hand, while the green orbit is not limited by the SZV, it does not escape from the system, due to stickiness phenomena in some unknown 2D surfaces. We plan to investigate these surfaces that prevent 3D chaotic orbits from escape, due to stickiness phenomena, in a future work.

Fig. 11. 3D PSS as shown in Fig. 9, together with the surface of zero velocity (SZV) for energy Ej = −4.0 seen from two different angles. The red chaotic orbit seems to be bounded by the SZV and cannot escape from the system, while the green chaotic orbit is not bounded by the SZV but is sticky to a hypersurface probably related to 2D unstable manifolds.

Next we wanted to determine the type and speed of diffusion of a sticky chaotic orbit finally escaping from the system. A common approach to quantifying chaotic diffusion and its parameters is through the measurement of the mean square displacement (MSD) (introduced by Zaslavsky 2002) of a random variable, $\mathbf{X}$, of an ensemble of n trajectories on the PSS. In our case, where the PSS is 4D, $\mathbf{X}=(\mathbf{y},\dot{\mathbf{y}},\mathbf{z},\dot{\mathbf{z}})$ (with x = 0 and $\dot{x}<0$). Therefore, the MSD, in the 4D phase space, is by definition

$$\begin{array}{cc}\hfill \mathit{MSD}& =⟨{(y-y_0)}^2⟩+⟨{(\dot{y}-{\dot{y}}_0)}^2⟩\hfill \\ \hfill & +⟨{(z-z_0)}^2⟩+⟨{(\dot{z}-{\dot{z}}_0)}^2⟩=2dD{t}^{\alpha },\hfill \end{array}$$(15)

where ⟨ ⋅ ⟩ denotes the mean value of n orbits in a very small grid around a specific point of the PSS ($y_0,\dot{y_0},z_0,\dot{z_0}$), d is the dimensionality of $\mathbf{X}$ (d = 4 in our case), D is the generalized diffusion coefficient, and α is the diffusion exponent that determines whether the process corresponds to normal diffusion (α = 1), to subdiffusion (0 ≤ α < 1 – slow diffusion), or to superdiffusion (1 < α – fast diffusion). The ballistic motion corresponds to a constant velocity with α = 2. In the case of normal diffusion, the only unknown parameter in Equation (15) is the diffusion coefficient, D, which can be determined numerically. In the case of anomalous diffusion, both D and α must be determined numerically.

It is more convenient to take the logarithm of both sides of Equation (15) and to fit a linear function. The slope of the fit then yields the diffusion exponent, α, whereas the vertical intercept gives the diffusion coefficient, D. In Fig. 12 the diffusion exponent, α, of Eq. (15) is plotted as a function of time for a sticky chaotic orbit (close and inside the blue orbit of Figs. 8, 9). We used the mean values of 20 orbits in a small grid around the central value ($y_0,\dot{y_0},z_0,\dot{z_0}$) in Eq. (15). For approximately the first 1400 iterations, which corresponds to approximately 14 Hubble times, α has a mean value around zero, denoting a regular behaviour, and then gets diffused outwards through a superdiffusing procedure. Finally, when the orbit escapes from the system, it follows a superballistic diffusion, following an Archimedean spiral (see Fig. 14, right) that seems to converge to a ballistic motion (with α→ 2 for time →∞). An Archimedean spiral is the locus corresponding to the locations over time of a point moving away from a fixed point with a constant speed along a line that rotates with constant angular velocity, so it is consistent with the ballistic motion.

Fig. 12. Diffusion exponent, α (Eq. (14)), as a function of time for a sticky chaotic orbit (close and inside the blue orbit of Figs. 8, 9). (a) For approximately the first 1400 iterations (approximately 14 Hubble times) the slope has a mean value around zero, denoting a regular behaviour due to stickiness phenomena, and then gets diffused outwards through a superdiffusing procedure. (b) After escaping from the system it presents a superdiffusing (or superballistic) behaviour and seems to converge to a ballistic motion (with α → 2 as t → ∞).

4.4. Apocentric sections

We wanted to investigate the morphological structures of the galaxy in the configuration space that are supported by the 3D sticky chaotic orbits that we studied in the previous section. In Figure 13, top, we plot the 2D PSS (y,ẏ) of planar orbits (in black), for the same energy level, E_j = −4.0, as in Figs. 6,8. We plot in colour the projections of the 4D PSS of 3D sticky chaotic orbits with initial conditions lying slightly above the unstable manifolds of each region (with z = 0.02). The red orbit lies above the unstable manifold of the planar PL₂ orbit. The red and magenta orbits stay located in these regions and behave like regular orbits for extremely long times, while the green and blue orbits escape from the system after a long time of integration, presenting the type of diffusion described in Fig. 12. This means that when an orbit escapes from the system it gets diffused very quickly (superdiffusion) outwards and then converges to a ballistic type of diffusion.

Fig. 13. Top: Projection (y,ẏ, with x = 0 and $\dot{x}<0$) of the 4D PSS of 3D sticky chaotic orbits having initially a small z component (z = 0.02) and lying above the planar unstable asymptotic manifolds of successive chaotic regions (in colour). The PSS of the planar orbits is in black. Bottom: Apocentric sections on the configuration space x − y of the 3D chaotic orbits corresponding to the same colours as in the top of this figure. They support the envelope of the bar as well as ring-like structures of R1 type that fill a small region around the bar. The planar unstable apocentric manifolds of orbits with initial conditions in the chaotic region of the 2D PSS outside the last KAM curve are in black. The energy level is Ej = −4.0.

In order to study the morphological structures of the galaxy that are supported by these 3D sticky chaotic orbits, we used the method of ‘apocentric manifolds’ (Voglis et al. 2006; Tsoutsis et al. 2008, 2009). In Figure 13, top, we plot the projection (y,ẏ, with x = 0 and $\dot{x}<0$) of the 4D PSS of 3D sticky chaotic orbits having initially a small z component (z = 0.02) and lying above the planar unstable asymptotic manifolds of successive chaotic regions (in colour), for energy level E_j = −4.0. The PSS of the planar orbits is in black. In Figure 13, bottom, we plot, in colour, only the apocentres of the projection of the 3D chaotic orbits described above, on the x − y configuration space (plane of rotation) using the same colours as in the figure above. We observe that the coloured apocentric sections of these 3D sticky chaotic orbits all support ring-like structures of R₁ type having a main axis perpendicular to the axis of the bar, for hundreds of Hubble times. The apocentric sections of the 3D sticky chaotic orbit above the planar PL₂ orbit (red) support the envelope of the bar as well as a ring structure. All of the different apocentric sections described above contribute to creating the thickness of the ring structure. The same is true for the stable periodic orbits lying among the chaotic regions, i.e. the 3D regular orbits lying above these periodic orbits (having initial conditions with a small perturbation on the z-axis) all support the ring structure. This is an outer ring, as it is placed close to the OLR, which is near 10 kpc.

In Figure 13, bottom, we plot in black the planar apocentric manifolds (the apocentres of all the orbits with initial conditions along the unstable manifold) of planar unstable periodic orbits lying outside the blue orbit. These form a path with spiral structures; however, this is not populated by 3D chaotic orbits, as orbits with initial conditions above these planar manifolds are rapidly diffused outwards (see blue apocentric sections) and escape from the system. Therefore, spiral structures are not supported in this case. This is due to the relatively small perturbation of the bar used in this model (only 12% in forces; see Fig. 1). We will study models with larger values of bar perturbation, M_b, in a future paper.

In Fig. 14 we plot the time evolution of the projection on the x − y plane of a 3D sticky chaotic orbit (blue in Fig. 13). The left panel corresponds to the first 10 000 iterations on the 4D PSS ($y,\dot{y},z,\dot{z}$), where the orbit behaves like a regular one, supports a ring-like structure, and exhibits subdiffusion with a diffusion exponent α = 0 (Fig. 12a). The panel in the middle corresponds to 20 000–35 000 iterations, where the orbit manifests its chaotic behaviour but is still bounded by the system, having a superdiffusive behaviour, while the panel on the right corresponds to 80 000–90 000 iterations, where the orbit escapes from the system through an Archimedean spiral and a ballistic diffusion (Fig. 12b).

Fig. 14. Time evolution of the projection on the x − y plane of a 3D sticky chaotic orbit (blue in Fig. 13). Left: First 10 000 iterations of the 4D PSS ($y,\dot{y},z,\dot{z}$). Middle: Iterations 20 000–35 000. Right: Iterations 80 000–90 000. The orbit finally escapes from the system through an Archimedean spiral.

The same phenomena happen in all energy levels above E_jL1, 2 (the energy of the L_1, 2 Lagrangian points). Figure 15 is the same as Fig. 13 for energy level E_j = −3.5. We again plot the PSS of the planar orbits in black, and in magenta we plot the (y,ẏ) projection of the 4D PSS of 3D sticky chaotic orbits initially having a small z component (z = 0.02) and lying above the planar unstable asymptotic manifolds of a chaotic region near the border of regular motion (in the top part of the figure). The planar PL₂ family is stable at this energy level. In the bottom part of the figure we plot in magenta the apocentric sections of the chaotic orbits plotted above, on the x − y plane. These support an outer ring-like structure of R₁ type that fills a small region around the bar. The planar unstable apocentric manifolds of orbits outside the last KAM curve are plotted in black. Here, again, we observe that the 3D sticky chaotic orbits do not support the spiral structures of the planar apocentric manifolds, but get quickly diffused outwards (with superdiffusion) once they escape from the ring structure.

Fig. 15. Top: The (y,ẏ, with x = 0 and $\dot{x}<0$) projection of the 4D PSS of 3D sticky chaotic orbits having a small z component (z = 0.02) initially and lying above the planar unstable asymptotic manifolds of a chaotic region near the border of regular motion (in magenta). The planar PL2 family is stable in this energy level. The PSS of the planar orbits is in black. Bottom: The apocentric sections of the chaotic orbits (in magenta) plotted above with the same colours, on the x − y plane. These support a ring-like structure of R1 type that fill a small region around the bar. The planar unstable manifold of orbits outside the last KAM curve is in black. The energy level is Ej = −3.5.

Figure 16 is the same as Fig. 15 for energy level E_j = −4.5. Here, we observe rings of R₁R₂′ type (a ring with its main axis perpendicular to the bar and a bigger pseudoring with its main axis parallel to the bar, which seems like a wrapped pair of spiral arms; see Athanassoula et al. 2009a,b) that are supported by the 3D sticky chaotic orbits. Here, the 3D sticky chaotic orbits do not support the spiral structures made from the planar apocentric manifolds, but get quickly diffused outwards once they escape from the ring structures.

Fig. 16. Same as Fig. 15, for energy level Ej = −4.5. Here, rings of R1R2′ type are supported by the 3D sticky chaotic orbits.

5. Conclusions and discussion

We study a Milky Way-like barred galactic model in order to understand how the internal dynamics define the morphological structures supported on the configuration space of the galaxy. Our model consists of a two-component axisymmetric disk and a bar potential with a rather small perturbation corresponding to a maximum of ≈12% in force perturbation, Q_f (see Fig. 1). We find the planar PL_1, 2 orbits that are the continuation of the Lagrangian equilibrium points L₁ and L₂ and exist for energies above E_L1, 2 (we also study their stability indices, b₁ and b₂) and we investigate the 4D PSS, giving priority to the 3D sticky chaotic orbits that support ring-like structures outside corotation for very long times compared to the Hubble time.

We study three different energy levels on the 4D PSS. We conclude that 3D sticky chaotic orbits with initial conditions just above the unstable manifolds of the planar PL₂ orbit, as well as other planar unstable periodic orbits, support ring-like structures of R₁ and R₁R₂ ′ type for extremely long times (hundreds and thousands of Hubble times) before escaping from the system.

We study the time evolution of the finite-time Lyapunov characteristic time (FT − LCN) as well as the speed and type of diffusion of these 3D sticky chaotic orbits. We conclude that they behave like regular orbits for extremely long times, compared to the Hubble time. Moreover, once these orbits escape from the system, they follow the same type of diffusion. They present subdiffusion during the time that they support the ring structures and then are quickly diffused (with superdiffusion) outwards. Once they escape from the system, they follow Archimedean spiral orbits, with a diffusion that tends to ballistic motion. Although there exist planar spiral arms on the configuration space made out of the apocentric manifolds of all the planar unstable periodic orbits outside corotation, these spiral paths are not populated by 3D chaotic orbits as they all get diffused very quickly outwards along these paths.

Buta and Combes 1996 claimed that, since the vast majority of ringed galaxies do not interact violently with other galaxies, it seems certain that rings are mainly a problem of internal dynamics. The present study gives an explanation of how the internal dynamics of a barred galactic model support ring-like structures. These results are also in agreement with Athanassoula et al. (2009a), who found from their simulations that less strong bars give rise to R₁ rings or pseudorings, while stronger bars drive R₂, R₁R₂, and spiral morphologies.

There are many observational studies that calculate the bar-induced perturbation strengths of barred spiral galaxies using the force perturbation, Q_f (see, for example, Block et al. 2004; Laurikainen et al. 2004; Buta et al. 2005, 2009). In these studies one can observe that barred galaxies with smaller force perturbation often present ring-like structures, while others with greater force perturbation present more open spiral arms.

Moreover, in Block et al. (2004), the authors verified a possible weak correlation between the strength of the bar and the strength of the spiral arms, in the sense that galaxies with the strongest bars tend to also have strong spirals with more open spiral arms (Athanassoula et al. 2009b). This is the topic of our future work (paper II), in which we are going to study the features supported by 3D sticky chaotic orbits in the case of stronger bar perturbation than tested here in our 3D galactic model. We will also examine the role of the pattern speed of the bar in supporting the morphological features of the galaxy.

Acknowledgments

This work was funded by the Sectoral Development Program (“OΠΣ” 5223471) of the Ministry of Education, Religious Affairs and Sports, through the National Development Program (NDP) 2021-25. It was conducted as part of project 200/1022, supported by the Research Committee of the Academy of Athens. We would like to thank Dr. Christos Efthymiopoulos for his help through our fruitful discussions on the subject of our paper.

References

Appendix A: Calculation of the stability indices of a periodic orbit

We consider a periodic orbit in a Hamiltonian system with three degrees of freedom, with initial conditions $y_0,z_0,{\dot{y}}_0,{\dot{z}}_0$ on the 4D Poincaré surface of section (PSS) x = 0 with $\dot{x}<0$ of the 6D phase space. The corresponding Poincaré map is T : ℜ⁴ → ℜ⁴. T connects the initial point $P(y_0,z_0,{\dot{y}}_0,{\dot{z}}_0)$ with the next point $P(y_1,z_1,{\dot{y}}_1,{\dot{z}}_1)$ on the corresponding PSS. This point is given from the equations:

$$\begin{array}{c}\hfill y_1=g_1(y_0,z_0,{\dot{y}}_0,{\dot{z}}_0)\\ \hfill z_1=g_2(y_0,z_0,{\dot{y}}_0,{\dot{z}}_0)\\ \hfill {\dot{y}}_1=g_3(y_0,z_0,{\dot{y}}_0,{\dot{z}}_0)\\ \hfill {\dot{z}}_1=g_4(y_0,z_0,{\dot{y}}_0,{\dot{z}}_0)\end{array}$$(A.1)

If we put a small perturbation ($\mathrm{\Delta }{y}_0,\mathrm{\Delta }{z}_0,\mathrm{\Delta }{\dot{y}}_0,\mathrm{\Delta }{\dot{z}}_0$) in the periodic orbit, the new point has coordinates ($y_0+\mathrm{\Delta }{y}_0,z_0+\mathrm{\Delta }{z}_0,{\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_0,{\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_0$). Then the Poincaré mapping of this point is:

$$\begin{array}{c}\hfill y_0+\mathrm{\Delta }{y}_1=g_1(y_0+\mathrm{\Delta }{y}_0,z_0+\mathrm{\Delta }{z}_0,{\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_0,{\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_0)\\ \hfill z_0+\mathrm{\Delta }{z}_1=g_2(y_0+\mathrm{\Delta }{y}_0,z_0+\mathrm{\Delta }{z}_0,{\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_0,{\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_0)\\ \hfill {\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_1=g_3(y_0+\mathrm{\Delta }{y}_0,z_0+\mathrm{\Delta }{z}_0,{\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_0,{\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_0)\\ \hfill {\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_1=g_4(y_0+\mathrm{\Delta }{y}_0,z_0+\mathrm{\Delta }{z}_0,{\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_0,{\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_0)\end{array}$$(A.2)

If we expand the previous equations (through Taylor series) we have:

$$\begin{array}{c}\hfill y_0+\mathrm{\Delta }{y}_1=y_0+\frac{\partial g_1}{\partial y}\mathrm{\Delta }{y}_0+\frac{\partial g_1}{\partial z}\mathrm{\Delta }{z}_0+\frac{\partial g_1}{\partial \dot{y}}\mathrm{\Delta }{\dot{y}}_0+\frac{\partial g_1}{\partial \dot{z}}\mathrm{\Delta }{\dot{z}}_0+...\\ \hfill z_0+\mathrm{\Delta }{z}_1=z_0+\frac{\partial g_2}{\partial y}\mathrm{\Delta }{y}_0+\frac{\partial g_2}{\partial z}\mathrm{\Delta }{z}_0+\frac{\partial g_2}{\partial \dot{y}}\mathrm{\Delta }{\dot{y}}_0+\frac{\partial g_2}{\partial \dot{z}}\mathrm{\Delta }{\dot{z}}_0+...\\ \hfill {\dot{y}}_0+\mathrm{\Delta }{\dot{y}}_1={\dot{y}}_0+\frac{\partial g_3}{\partial y}\mathrm{\Delta }{y}_0+\frac{\partial g_3}{\partial z}\mathrm{\Delta }{z}_0+\frac{\partial g_3}{\partial \dot{y}}\mathrm{\Delta }{\dot{y}}_0+\frac{\partial g_3}{\partial \dot{z}}\mathrm{\Delta }{\dot{z}}_0+...\\ \hfill {\dot{z}}_0+\mathrm{\Delta }{\dot{z}}_1={\dot{z}}_0+\frac{\partial g_4}{\partial y}\mathrm{\Delta }{y}_0+\frac{\partial g_4}{\partial z}\mathrm{\Delta }{z}_0+\frac{\partial g_4}{\partial \dot{y}}\mathrm{\Delta }{\dot{y}}_0+\frac{\partial g_4}{\partial \dot{z}}\mathrm{\Delta }{\dot{z}}_0+...\end{array}$$(A.3)

If we keep only the terms of the first order of the Taylor series in eq. (A.3) we have a linear system for the Poincaré mapping of the perturbations:

$\left[\begin{array}{c}\mathrm{\Delta }{y}_1\\ \mathrm{\Delta }{z}_1\\ \mathrm{\Delta }{\dot{y}}_1\\ \mathrm{\Delta }{\dot{z}}_1\end{array}\right]=\left[\begin{array}{c}{a}_{11}{a}_{12}{a}_{13}{a}_{14}\\ a_{21}{a}_{22}{a}_{23}{a}_{24}\\ a_{31}{a}_{32}{a}_{33}{a}_{34}\\ a_{41}{a}_{42}{a}_{43}{a}_{44}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{y}_0\\ \mathrm{\Delta }{z}_0\\ \mathrm{\Delta }{\dot{y}}_0\\ \mathrm{\Delta }{\dot{z}}_0\end{array}\right]$

where

$$\begin{array}{c}\hfill a_{i1}=\frac{\partial g_i}{\partial y},a_{i2}=\frac{\partial g_i}{\partial z},,a_{i3}=\frac{\partial g_i}{\partial \dot{y}},a_{i4}=\frac{\partial g_i}{\partial \dot{z}}\end{array}$$(A.4)

with i = 1, ..., 4. This is equivalent to the expression:

$$\begin{array}{c}\hfill \mathit{\xi }=A{\mathit{\xi }}_0\end{array}$$(A.5)

where ξ₀ is the initial deviation vector of the periodic orbit, ξ is the deviation vector after the first Poincaré mapping and A is a 4 × 4 matrix called the monodromy matrix of T.

In order to determine the kind of stability of the periodic orbit (stable or unstable) we find the 4 eigenvalues of the monodromy matrix A through the equation:

$$\begin{array}{c}\hfill det(A-\lambda I)=0\end{array}$$(A.6)

where I is the 4 × 4 unitary matrix. The 4 eigenvalues are found from the equation:

$$\begin{array}{c}\hfill {\lambda }^4+a{\lambda }^3+\beta {\lambda }^2+a\lambda +\delta =0\end{array}$$(A.7)

The coefficients of the polynomial are (see Contopoulos 2002, p.286):

$$\begin{array}{c}\hfill a=-(a_{11}+a_{22}+a_{33}+a_{44})=-tr(A)\end{array}$$(A.8)

$$\beta =\left|\begin{array}{c}{a}_{11}{a}_{12}\\ a_{21}{a}_{22}\end{array}\right|+\left|\begin{array}{c}{a}_{11}{a}_{13}\\ a_{31}{a}_{33}\end{array}\right|+\left|\begin{array}{c}{a}_{11}{a}_{14}\\ a_{41}{a}_{44}\end{array}\right|+\left|\begin{array}{c}{a}_{22}{a}_{23}\\ a_{32}{a}_{33}\end{array}\right|+\left|\begin{array}{c}{a}_{22}{a}_{24}\\ a_{42}{a}_{44}\end{array}\right|+\left|\begin{array}{c}{a}_{33}{a}_{34}\\ a_{43}{a}_{44}\end{array}\right|$$

$$\begin{array}{c}\hfill \delta =det(A)=1\end{array}$$(A.9)

The last equation (A.9) is a consequence of the preservation of the volumes in the 4-D phase space (Hadjidemetriou 1975).

A.1. The numerical computation of the elements of the monodromy matrix of a periodic orbit:

We perturb (with a perurbation of the order Δy₀ = ϵ = 10⁻⁸ or 10⁻⁹) the initial condition y₀ of the periodic orbit (y₀, z₀, ẏ₀, $\dot{z_0}$) and we integrate using a Runge Kutta method in order to find the next iteration on the Poincaré surface of section: (y₀ + Δy₀, z₀, ẏ₀, $\dot{z_0}$) → (y₁,z₁, $\dot{y_1}$,$\dot{z_1}$), with x = 0 and ẋ < 0. This intersection gives us the first four elements of the monodromy matrix a_ij, i = 1, ..., 4 and j = 1 through the following equations:

$$\begin{array}{cc}\hfill a_{11}& =\frac{y_1-y_0}{ϵ}\hfill \\ \hfill a_{21}& =\frac{z_1-z_0}{ϵ}\hfill \\ \hfill a_{31}& =\frac{{\dot{y}}_1-{\dot{y}}_0}{ϵ}\hfill \\ \hfill a_{41}& =\frac{{\dot{z}}_1-{\dot{z}}_0}{ϵ}\hfill \end{array}$$(A.10)

Consequently, we perturb the initial condition z₀ of the periodic orbit with a perturbation Δz₀ and we compute the next intersection on the Poincaré surface of section: (y₀, z₀ + Δ_z0, ẏ₀, $\dot{z_0}$) → (y₂,z₂, $\dot{y_2}$,$\dot{z_2}$) with x = 0 and $\dot{x}<0$. This intersection gives us the next four elements of the monodromy matrix, a_ij, i = 1, ..., 4 and j = 2. Similarly, we perturb the initial conditions ${\dot{y}}_0$ and ${\dot{z}}_0$ and we find the elements a_ij, i = 1, ..., 4 with j = 3 and j = 4 respectively.

A.2. The stability indices b₁ and b₂

The stability indices of the periodic orbit are defined as:

$$\begin{array}{c}\hfill b_{1,2}=\frac{a\pm \sqrt{\mathrm{\Delta }}}{2}\end{array}$$(A.11)

where:

$$\begin{array}{c}\hfill \mathrm{\Delta }=a^2-4(\beta -2)\end{array}$$(A.12)

The eigenvalues of the periodic orbit are computed through the stability indices (see Broucke 1969; Hadjidemetriou 1975):

$$\begin{array}{c}\hfill {\lambda }_1,\frac{1}{{\lambda }_1}=\frac{-b_1\pm \sqrt{b_1^2-4}}{2}\\ \hfill {\lambda }_2,\frac{1}{{\lambda }_2}=\frac{-b_2\pm \sqrt{b_2^2-4}}{2}\end{array}$$(A.13)

The periodic orbit is:

a) Stable if Δ > 0 and |b_1, 2|< 2.

b) Simply unstable if Δ > 0 and |b₁|< 2 and |b₂|> 2 or |b₁|> 2 and |b₂|< 2.

c) Double unstable if Δ > 0 and |b₁|> 2 and |b₂|> 2.

d) Complex unstable if Δ < 0.

Appendix B: Calculation of the Lyapunov Characteristic number (LCN)

In order to determine the chaotic nature of orbits we use the numerical computation of the Lyapunov Characteristic Number (LCN) χ:

$$\begin{array}{c}\hfill \chi =\underset{t\to \infty }{lim}\frac{1}{t}ln\frac{w(t)}{w(0)}\end{array}$$(B.1)

where w(0) is an initial deviation vector of the orbit and w(t) is the deviation vector after a time t. For a regular orbit χ = 0, while for a chaotic one χ > 0, implying exponential divergence of nearby orbits. The computation of the LCN has been used extensively as a chaos indicator after the introduction of numerical algorithms for the determination of its value at late 70’s [Benettin et al. 1976; Benettin & Strelcyn 1978; Contopoulos et al. 1978].

In order to calculate the LCNχ of an orbit of the 3D Hamiltonian (8) we use a system of 6 variational equations (Skokos 2010):

$$\begin{array}{c}\hfill \dot{\mathit{w}}=[{\mathit{J}}_6·{\mathit{D}}^{\mathbf{2}}\mathit{H}(x(t))]·w\end{array}$$(B.2)

where $\dot{\mathit{w}}$ is the time evolution of an initial deviation vector w(0) = (δx(0),δy(0),δz(0),δp_x(0),δp_y(0),δp_z(0)), D²H(x(t)) is the Hessian matrix of Hamiltonian (8) calculated on the reference orbit x(t) = (x, y, z, p_x, p_y, p_z) and J₆= $\left[\begin{array}{c}{0}_3{I}_3\\ -I_3{0}_3\end{array}\right]$.

The system of the 6 variational equations, following eq. (B.2), for a specific orbit is the following:

$$\begin{array}{c}\hfill \dot{\delta x}=\frac{{\partial }^{2}H}{\partial p_x\partial x}\delta x+\frac{{\partial }^{2}H}{\partial p_x\partial y}\delta y+\frac{{\partial }^{2}H}{\partial p_x\partial z}\delta z+\frac{{\partial }^{2}H}{\partial p_x^2}\delta p_x+\frac{{\partial }^{2}H}{\partial p_x\partial p_y}\delta p_y+\frac{{\partial }^{2}H}{\partial p_x\partial p_z}\delta p_z\\ \hfill \dot{\delta y}=\frac{{\partial }^{2}H}{\partial p_y\partial x}\delta x+\frac{{\partial }^{2}H}{\partial p_y\partial y}\delta y+\frac{{\partial }^{2}H}{\partial p_y\partial z}\delta z+\frac{{\partial }^{2}H}{\partial p_y\partial p_x}\delta p_x+\frac{{\partial }^{2}H}{\partial p_y^2}\delta p_y+\frac{{\partial }^{2}H}{\partial p_y\partial p_z}\delta p_z\\ \hfill \dot{\delta z}=\frac{{\partial }^{2}H}{\partial p_z\partial x}\delta x+\frac{{\partial }^{2}H}{\partial p_z\partial y}\delta y+\frac{{\partial }^{2}H}{\partial p_z\partial z}\delta z+\frac{{\partial }^{2}H}{\partial p_z\partial p_x}\delta p_x+\frac{{\partial }^{2}H}{\partial p_z\partial p_y}\delta p_y+\frac{{\partial }^{2}H}{\partial p_z^2}\delta p_z\\ \hfill \dot{\delta p_x}=-\frac{{\partial }^{2}H}{\partial x^2}\delta x-\frac{{\partial }^{2}H}{\partial x\partial y}\delta y-\frac{{\partial }^{2}H}{\partial x\partial z}\delta z-\frac{{\partial }^{2}H}{\partial x\partial p_x}\delta p_x-\frac{{\partial }^{2}H}{\partial x\partial p_y}\delta p_y-\frac{{\partial }^{2}H}{\partial x\partial p_z}\delta p_z\\ \hfill \dot{\delta p_y}=-\frac{{\partial }^{2}H}{\partial y\partial x}\delta x-\frac{{\partial }^{2}H}{\partial y^2}\delta y-\frac{{\partial }^{2}H}{\partial y\partial z}\delta z-\frac{{\partial }^{2}H}{\partial y\partial p_x}\delta p_x-\frac{{\partial }^{2}H}{\partial y\partial p_y}\delta p_y-\frac{{\partial }^{2}H}{\partial y\partial p_z}\delta p_z\\ \hfill \dot{\delta p_z}=-\frac{{\partial }^{2}H}{\partial z\partial x}\delta x-\frac{{\partial }^{2}H}{\partial z\partial y}\delta y-\frac{{\partial }^{2}H}{\partial z^2}\delta z-\frac{{\partial }^{2}H}{\partial z\partial p_x}\delta p_x-\frac{{\partial }^{2}H}{\partial z\partial p_y}\delta p_y-\frac{{\partial }^{2}H}{\partial z\partial p_z}\delta p_z\end{array}$$(B.3)

When integrating an orbit we use simultaneously the 6 equations of motion (9) and the 6 variational equations (B.3) for finding the time evolution of the deviation vector w(t) = (δx²(t)+δy²(t)+δz²(t)+δp_x²(t)+δp_y²(t)+δp_z²(t))^1/2, of an initial deviation vector w(0) and we calculate the Lyapunov characteristic number form eq. (B.1). The LCN is calculated on the 4D PSS $y,\dot{y},z,\dot{z}$ for x = 0 and $\dot{x}<0$. We use initial values of the deviations: δx(0) = δy(0) = δz(0) = δp_x(0) = δp_y(0) = δp_z(0) = 1 and then we perform a normalization of the deviation vector in each step of the integration in order to avoid numerical overflow.

The inverse of the LCN, which is called Lyapunov time: t_L = 1/χ gives an estimate of the time needed for a dynamical system to become chaotic and in practice measures the time needed for nearby orbits of the system to acquire exponential divergence.

All Figures

	Fig. 1. Q-force perturbation of the bar, Qf, as a function of the radius, r, of the bar. The maximum Qf is approximately 12%. The vertical line gives the position of the Lagrangian points L1, 2 that are located shortly after the end of the bar.
In the text

Fig. 2. (a) Stability indices, b1, b2, of the planar PL1, 2 orbit (branch 1 and branch 2) as a function of the Jacobi constant (or energy), En, for bar mass Mb = 5. The vertical stability index, b1, (black) corresponds to the stability of planar periodic orbits with respect to perturbations perpendicular to the orbital plane, while the planar stability index, b2, (red) corresponds to the stability with respect to perturbations on the orbital plane. The three vertical lines correspond to the energy levels at which we investigated the phase space. (b) Characteristics of the two planar branches of the PL1 family that are joined and disappear at the bifurcation point. (c) Branch 1 and (d) branch 2 planar PL1 orbits plotted with the same colour as the corresponding energies (vertical lines) in Fig. 2a.

In the text

Fig. 3. 4D PSS ($y,\dot{y},z,\dot{z}$) for x = 0 and $\dot{x}<0$, for energy level Ej = −4.5. (a) For smaller radii. (b) For larger radii. The colour gradient allows the visualization of the fourth dimension, ż, in rotational tori of regular orbits according to the colour bar on the left (minimum ż corresponds to 0 and maximum ż to 1). The rotational torus of the regular ring-like orbit is also shown. The 2D PSS of the planar orbits is plotted in black and the unstable manifolds of the planar unstable periodic orbits in the chaotic region are in red. 3D sticky chaotic orbits (green) escape from the system after very long times. (a′) and (b′) are the 2D projections (y,ẏ) of the 4D PSS of the same orbits as in (a) and (b).

In the text

	Fig. 4. Same as in Fig. 3, but for energies (a), (a′) Ej = −4.0 and (b), (b′) Ej = −3.5.
In the text

	Fig. 5. Top: 3D regular orbit that supports an R1 ring-like structure of the galaxy, labelled ‘Ring Orbit’ (red), together with the planar unstable periodic orbit (black). Bottom: 2D PSS (y,ẏ) for x = 0 and $\dot{x}<0$. The projection of the 3D ring-like orbit is plotted in red. The energy level is Ej = −4.5.
In the text

	Fig. 6. 2D PSS (y,ẏ) for x = 0 and $\dot{x}<0$ of planar orbits (with $z,\dot{z}=0$) for energy Ej = −4.0. The coloured islands of stability correspond to orbits shown in Fig. 7.
In the text

	Fig. 7. Planar orbits corresponding to the islands of stability of the same colours as in Fig. 6. These all support ring-like structures.
In the text

Fig. 8. 2D PSS (y,ẏ, for x = 0 and $\dot{x}<0$) of planar orbits for energy level Ej = −4.0 (black). The (y,ẏ) projections of the 4D PSS of the first 20 000 iterations of three 3D sticky chaotic orbits having a small third dimension z = 0.02 above the plane of rotation are plotted in colour. The red one behaves like a regular orbit and the green one like a weakly chaotic orbit for very long times. The blue one is diffused rather quickly outwards.

In the text

Fig. 9. Top: 3D projection of the 4D PSS ($y,\dot{y},z$) of the 3D chaotic orbits plotted in colour in Fig. 8. Top left: Orbits for the first 20 000 iterations. The red and green orbits behave as ordered ones for very long times, while the blue one is diffused out of the system more rapidly. Top middle: Orbits for iterations between 20 000 and 50 000. Top right: Orbits for the next 100 000 iterations. The time evolutions of the red and green sticky chaotic orbits present a very slow diffusion inwards. Bottom: Corresponding 2D projections of the 4D PSS for the same time periods. The 2D PSS of the planar orbits is shown in black, with the projections of the sticky 3D chaotic orbits presented above shown in colour.

In the text

Fig. 10. Top: 4D PSS for energy level Ej = −4.0. The same as in Fig. 9, but here in orange or yellow (the gradient of the colour indicates the ż component) showing a sticky chaotic orbit that acts as an approximate barrier between the orbits that get diffused inwards (green) and those that get diffused outwards (blue). Bottom: Time evolution of FT-LCN for the two sticky chaotic orbits of Fig. 9. A chaotic orbit (blue) diffused faster outwards and reaching very large distances, and a sticky chaotic orbit (green) that behaves like a regular one for very long times and finally converges to a positive LCN.

In the text

	Fig. 11. 3D PSS as shown in Fig. 9, together with the surface of zero velocity (SZV) for energy Ej = −4.0 seen from two different angles. The red chaotic orbit seems to be bounded by the SZV and cannot escape from the system, while the green chaotic orbit is not bounded by the SZV but is sticky to a hypersurface probably related to 2D unstable manifolds.
In the text

Fig. 12. Diffusion exponent, α (Eq. (14)), as a function of time for a sticky chaotic orbit (close and inside the blue orbit of Figs. 8, 9). (a) For approximately the first 1400 iterations (approximately 14 Hubble times) the slope has a mean value around zero, denoting a regular behaviour due to stickiness phenomena, and then gets diffused outwards through a superdiffusing procedure. (b) After escaping from the system it presents a superdiffusing (or superballistic) behaviour and seems to converge to a ballistic motion (with α → 2 as t → ∞).

In the text

Fig. 13. Top: Projection (y,ẏ, with x = 0 and $\dot{x}<0$) of the 4D PSS of 3D sticky chaotic orbits having initially a small z component (z = 0.02) and lying above the planar unstable asymptotic manifolds of successive chaotic regions (in colour). The PSS of the planar orbits is in black. Bottom: Apocentric sections on the configuration space x − y of the 3D chaotic orbits corresponding to the same colours as in the top of this figure. They support the envelope of the bar as well as ring-like structures of R1 type that fill a small region around the bar. The planar unstable apocentric manifolds of orbits with initial conditions in the chaotic region of the 2D PSS outside the last KAM curve are in black. The energy level is Ej = −4.0.

In the text

	Fig. 14. Time evolution of the projection on the x − y plane of a 3D sticky chaotic orbit (blue in Fig. 13). Left: First 10 000 iterations of the 4D PSS ($y,\dot{y},z,\dot{z}$). Middle: Iterations 20 000–35 000. Right: Iterations 80 000–90 000. The orbit finally escapes from the system through an Archimedean spiral.
In the text

Fig. 15. Top: The (y,ẏ, with x = 0 and $\dot{x}<0$) projection of the 4D PSS of 3D sticky chaotic orbits having a small z component (z = 0.02) initially and lying above the planar unstable asymptotic manifolds of a chaotic region near the border of regular motion (in magenta). The planar PL2 family is stable in this energy level. The PSS of the planar orbits is in black. Bottom: The apocentric sections of the chaotic orbits (in magenta) plotted above with the same colours, on the x − y plane. These support a ring-like structure of R1 type that fill a small region around the bar. The planar unstable manifold of orbits outside the last KAM curve is in black. The energy level is Ej = −3.5.

In the text

	Fig. 16. Same as Fig. 15, for energy level Ej = −4.5. Here, rings of R1R2′ type are supported by the 3D sticky chaotic orbits.
In the text