{'fraction_non_alphanumeric': 0.038798300510995004, 'fraction_numerical': 0.027372681862548084, 'mean_word_length': 4.162196206283343, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Instabilities and stickiness in a 3D rotating galactic potential

We study the dynamics in the neighborhood of simple and double unstable periodic orbits in a rotating 3D autonomous Hamiltonian system of galactic type. In order to visualize the four dimensional spaces of section we use the method of color and rotation. We investigate the structure of the invariant manifolds that we found in the neighborhood of simple and double unstable periodic orbits in the 4D spaces of section. We consider orbits in the neighborhood of the families x1v2, belonging to the x1 tree, and the z-axis (the rotational axis of our system). Close to the transition points from stability to simple instability, in the neighborhood of the bifurcated simple unstable x1v2 periodic orbits we encounter the phenomenon of stickiness as the asymptotic curves of the unstable manifold surround regions of the phase space occupied by rotational tori existing in the region. For larger energies, away from the bifurcating point, the consequents of the chaotic orbits form clouds of points with mixing of color in their 4D representations. In the case of double instability, close to x1v2 orbits, we find clouds of points in the four dimensional spaces of section. However, in some cases of double unstable periodic orbits belonging to the z-axis family we can visualize the associated unstable eigensurface. Chaotic orbits close to the periodic orbit remain sticky to this surface for long times (of the order of a Hubble time or more). Among the orbits we studied we found those close to the double unstable orbits of the x1v2 family having the largest diffusion speed.


Introduction
The aim of this paper is to study the dynamics in the neighborhood of unstable periodic orbits in a rotating three dimensional (hereafter 3D) autonomous Hamiltonian system of galactic type. Instability of periodic orbits in such a dynamical system is of three kinds: simple, double and complex instability. The case of complex instability was studied in a recent paper by Katsanikas et al ).
Here we study the two other cases. Furthermore, we examine the stickiness effects in our 3D system. This is the first time that the stickiness phenomenon is explicitly investigated in association with these kinds of instability (simple and double) in a 3D autonomous Hamiltonian system. For this purpose we consider the intersections of the orbits by a space of section. In a 3D autonomous Hamiltonian system the space of section has 4 dimensions. So, for the visualization of the 4D spaces of section it is convenient to use the method of color and rotation (Patsis & Zachilas 1994). Rotation is used for understanding the geometry of a 3D projection and color is used to visualize the 4th dimension.
3D rotating autonomous Hamiltonian systems are of special importance in Galactic Dynamics. They represent in a realistic way rotating galaxies, in which rotation plays a key role in the dynamics on the equatorial plane, as well as in the third dimension. The most important family for the vertical morphology of rotating disk galaxies is a family called x1v1 by . The x1v1 family is introduced at the transition from stability to simple instability (S → U ) as the stability of the central, planar family x1 undergoes a S → U → S transition within a small Jacobi integral interval at the vertical 2:1 resonance region (for a definition of resonances in rotating galactic potentials see Contopoulos pages 379 and 422). At the U → S transition, another family of periodic orbits (hereafter p.o.), x1v2, is bifurcated from x1 as simple unstable . For larger Jacobi constants x1v2 becomes double unstable. Thus, it possesses both kinds of instabilities we want to study. In this paper we focus on x1v2 and we examine the structure of the phase space in its neighborhood, when it is simple or double unstable. However, in order to investigate in depth the dynamics in the neighborhood of double unstable periodic orbits we will also study a case referring to z-axis orbits.
The importance of the dynamical phenomena at the 2:1 vertical resonance is due to the role of the associated families in shaping the vertical morphology of the 3D disks. They form to a large extent the peanut-shaped structure, which is characteristic of many disk galaxies observed edge-on . In addition, chaotic motion in the neighborhood of unstable p.o.'s has been proven to be significant for the dynamics of rotating galactic disks on their equatorial planes (for a review see e.g. Contopoulos 2009). In that respect the knowledge of the structure of phase space in the neighborhood of 3-dimensional p.o. in galactic systems is the first step for the investigation of a possible contribution of chaotic orbits to the observed vertical profiles of this type of galaxies.
The space of section that we use in this paper for investigating the orbital behavior close to x1v2 orbits is defined by y = 0 and we consider intersections withẏ > 0. Hereafter by "intersections" we will refer to intersections, which fulfill these properties, until we start discussing the dynamics close to z-axis. Then we will change space of section. Every consequent in our cross section is determined by 4 Cartesian coordinates, e.g. (x,ẋ, z,ż). In the method of color and rotation we project the consequents in a 3D subspace and observe them from different viewpoints. In the text the viewpoints are given in spherical coordinates. The unit for the distance d of the consequents of the surface of section from the observer corresponds to the longest side of the bounding box of the figure (more details can be found in Katsanikas and Patsis 2010). In practice we rotate the 3D figure in the computer screen in order to find the most appropriate viewpoint for the presentation on the 2D surface of the pages of a paper. We select the viewpoint that represents best the morphology of the structures formed by the consequents in the 3D projections of the 4D space of section in order to facilitate the reader to follow our description. Then we color every point according the value of its 4th coordinate. For this reason we use the "Mathematica" package and we map the interval of the values of the fourth coordinate into [0, 1].
• A well defined structure in all 3D projections and a smooth color variation along them is evidence for regular motion.

Stability of periodic orbits
The linear stability of the periodic orbits is studied, following the analysis by Broucke (Broucke 1969) in the case of the elliptic restricted three body problem and the method used by Hadjidemetriou (Hadjidemetriou 1975) in the three body problem. The method has been also described in a recent earlier paper of us in the series of papers studying the dynamics in the neighborhood of p.o. in 3D rotating galactic potentials (Katsanikas & Patsis 2010). The characteristic equation in a 3D autonomous hamiltonian system is of the form λ 4 + aλ 3 + βλ 2 + aλ + 1 = 0. and has 4 eigenvalues λ 1 , 1/λ 1 , λ 2 , 1/λ 2 , where λ i , 1 and ∆ = a 2 − 4(β − 2). Thus, the stability of the p.o. depends on the two stability indices b 1 and b 2 and on the quantity ∆.
Then we have three kinds of instability: (1) simple instability (U ): ∆ > 0 and |b 1 | > 2, |b 2 | < 2 or |b 1 | < 2, |b 2 | > 2, (two eigenvalues on the unit circle and two on the real axis) (2) double instability (DU ): ∆ > 0 and |b 1 | > 2 and |b 2 | > 2 (four eigenvalues on the real axis) (3) complex instability (∆) (all four complex eigenvalues are off the unit circle in the complex plane): ∆ < 0 (for details see Contopoulos 2002, p.286). Hereafter we will indicate the periodic orbits as S, U , DU and ∆ according to their stability. A generalization of the three kinds of instability in Hamiltonian systems of N > 3 degrees of freedom can be found in Skokos (Skokos 2001).
In section 2 we present our Hamiltonian system. In section 3 we discuss the stability of a family of p.o. we call x1v2. In section 4 we study the orbital behavior close to simple unstable orbits of x1v2, in section 5 we investigate the dynamical behavior in the neighborhood of double unstable orbits of x1v2, as well as of the z-axis families, in section 6 we compare the diffusion of the chaotic orbits away from the p.o. and finally in section 7 we discuss and enumerate our conclusions.

The Hamiltonian System
The Hamiltonian of our system in Cartesian coordinates is: where Φ(x, y, z) is the potential: This potential in its axisymmetric form (q a = 1, q b = 1) can be considered as an approximation of the potential of the Milky Way (Miyamoto & Nagai 1975).
In our units, the distance R=1 corresponds to 1 kpc. The velocity unit corresponds to 209.64 km/sec. The Jacobi constant Ej = 1 corresponds to 43950 (km/sec) 2 and we will call it "energy". For the parameters we have used the following values: a 1 = 0 kpc, b 1 = 0.495 kpc, M 1 = 2.05 × 10 10 M ⊙ , a 2 = 7.258 kpc, b 2 = 0.520 kpc, M 2 = 25.47 × 10 10 M ⊙ , q a = 1.2, q b = 0.9. For the angular velocity we used the value Ω b = 0.045 in our units. The parameters q a , q b determine the geometry of the disks, while a, b are scaling factors. The specific values of the geometrical parameters do not have a particular physical meaning. The dynamical phenomena we study in our paper are typical for rotating triaxial systems. The Ω b value we use here is different from the one we used in (Katsanikas & Patsis 2010) in order to have larger regions in the parameters' space with non-escape orbits close to double unstable periodic orbits (Section 5).

Stability Diagram
In order to study the evolution of the stability of a family of periodic orbits as a parameter of our system varies we use the stability diagram (Contopoulos 2002 p.287). The stability diagram gives the values of the coefficients b 1 and b 2 for successive periodic orbits in a family by changing only one parameter of the system (in our case E j ), while keeping the other parameters constant. In Fig. 1 we observe in such a diagram how the family of periodic orbits x1v2  bifurcates from the 2D central family x1 on the plane z = 0 at the U → S transition that happens in the vertical 2:1 resonance region (Contopoulos 2002 p.390). This bifurcation occurs for E j ≃ −4.7. We observe that the stability parameters b 1 , b 2 of the family x1v2 decrease as E j increases. The value of b 2 is initially about 1.95 at E j = −4.7 (left side of diagram). When b 2 becomes −2 we have a transition of x1v2 from simple to double instability (U → DU ) (for E j ≃ −3.41) and the x1v2 periodic orbits become double unstable. The E j interval of x1v2, where it is simple unstable (−4.7 < E j < −3.41) and the E j interval of x1v2 where it is double unstable (−3.41 < E j < −2.3), will be called hereafter regions U and DU respectively (Fig. 1).

Simple Instability
First we study the orbital behavior of the x1v2 family in the region U. In the following applications we use the set of projections (x,ẋ, z) and (ẋ, z,ż). However, we find qualitatively similar results by using the (x,ẋ,ż) and (x, z,ż) projections respectively. The use of the two projections gives us all needed information to understand the structure of the distribution of the consequents in the 4D space by means of the colorrotation method, as we will see below.
The initial conditions to be studied are chosen in two different ways. First we investigate the dynamical behavior of the orbits in the neighborhood of the simple unstable periodic orbit x1v2 at E j = −4.66 by varying the values and the direction (∆x, ∆ẋ, ∆z or ∆ż) of the perturbation. Secondly, we consider a constant perturbation ∆x = 10 −4 and we apply this perturbation to the initial conditions of several periodic orbits of the x1v2 family in the region U for different E j values. We start by keeping E j constant and investigate as a sample case the dynamical behavior very close to the point where we have the U → S transition of the x1 family and the bifurcation of the family x1v2 as U .

Orbital behavior close to the bifurcating point -A sample case
The energy E j = −4.66 is indicated with an arrow in Fig. 1. The periodic orbit has initial conditions (x 0 ,ẋ 0 , z 0 ,ż 0 ) = (0.037941612, 0, 0, 0.42902353). In order to study the orbital behavior in its neighborhood we added a small perturbation in theż-direction equal to ∆ż = −6 × 10 −6 . We call this perturbed orbit "u". The first 2580 consequents of orbit "u" define a ribbon of the form of a double loop in the (x,ẋ, z) subspace (Fig. 2). In the bottom of the figure we give also from left to right the 2D projections (x,ẋ), (x, z) and (ẋ, z) in order to help the reader understand the topology of the structure. The two lobes of the double loop are inclined with respect to each other and form an angle around 50 o in the (x, z) projection (middle panel in the lower part of Fig. 2). The (x, z) projection clearly also shows that the thickness of the lobes is very small (but not zero) and this is the reason we describe them as ribbons. Essentially the branches of the lobes are thin tubes. The 4D representation of this structure, colored according to the normalizedż values and using (x,ẋ, z) for the 3D spatial projection is given in Fig. 2 (top). For the 3D projection we use here a point of view defined by the angles (θ, φ) = (30 o , 108 o ). We observe that this ribbon has a smooth color variation from red to violet, in agreement with structures found in the neighborhood of simple periodic (U ) orbits by Patsis and Zachilas (1994). The structure of the phase space close to simple unstable periodic orbits has been studied for the first time by Magnenat (1982). However, Magnenat used 2D projections and these did not allow the full understanding of the 4D structure of the double loop objects. In the present case, our 4D representation allows us to understand the details of this structure.
By inspection of Fig. 2 we infer that the color variation along the two branches is such that at the intersection region we have the meeting of two different shades. This means that in this region the points have different values in the 4th dimension and there is no real intersection in the 4D space of section. Consequently the ribbon-like double loop structure has not an intersection in 4D. We observe that the colors at the region of the intersection are not only different, but that they are at the two extremes of our color scale, i.e. violet and red (see right side of Fig. 2). The colors in this region are clustered around two values (ż = ±0.429) very close to theż initial condition of the periodic orbitż = 0.42902353. Figure 2 reflects the symmetry of the galactic potential we study. Due to this symmetry, the families of p.o. that bifurcate at the vertical resonances from the planar x1 family come in pairs with mirror morphology with respect to the equatorial plane of the system in the configuration space, i.e. x1v1 and x1v1 ′ , x1v2 and x1v2 ′ etc. The p.o. of x1v2 and x1v2 ′ at every E j differ only in the sign of theirż 0 initial conditions. In the 3D (x,ẋ, z) projection we use in Fig. 2, the initial conditions of x1v2 and x1v2 ′ are identical and the p.o. is in the center of the intersection region of the two branches of the double loop ribbon (indicated with an arrow). The clustering of the colors around the two values (ż 0 ≈ ±0.429) in the central region is due to the coexistence of x1v2 and x1v2 ′ . The consequents drift away from the central region of the double loop by deviating from the two (±ż 0 ) initial conditions with different colors. The first 2580 consequents are along the double loop structure, thus they indicate a regular rather than chaotic motion.
We reach the same conclusion by using other 3D projections for the spatial coordinates of the consequents. In Fig. 3 the first 2580 consequents of the orbit "u" in the 3D projection (ẋ, z,ż) form again a ribbon-like surface, however this time with a single loop morphology. Color is determined by the x value of the consequents and we observe a smooth color variation along the ribbon. In this case the two representatives of the p.o x1v2 and x1v2 ′ are at two different positions in the (ẋ, z,ż) space, since they have differentż 0 initial conditions. Their positions are indicated with arrows in Fig. 3. The regions close to both of them are colored red (same x values) and departing from them, red becomes orange, then yellow, green, light blue, blue and finally violet. The ribbon is a 4D object. From both Figs. 2 and 3 we conclude that the first 2580 consequents form a well-defined ribbon-like surface that has no self-intersections in 4D.
A similar analysis can be done by means of the two other 3D projections (x,ẋ,ż) and (x, z,ż), which give 4D structures morphologically similar to (x,ẋ, z) and (ẋ, z,ż) respectively. It is evident that the orbit "u" in the neighborhood of a U x1v2 p.o. for 2580 consequents is confined in a particular region of the space of section with dimensions (∆x × ∆ẋ × ∆z × ∆ż) < (0.00035 × 0.0005 × 0.04 × 1). As we will see below, for larger integration times the orbit diffuses to a larger region of the phase space. The structures we find in the phase space indicate a sticky behavior and this will be further analyzed in section 4.4.
We note that although in general we have different colors in the intersections of double loop type structures, like the one in Fig. 2a, in exceptional cases we can have also real intersections in the 4D space. Such a case is encountered in the neighborhood of U z-axis orbits in 4D spaces of section defined by z = 0. The z-axis orbits lie entirely on the rotational, z axis, of the system. They have initial conditions  (x,ẋ, y,ẏ) = (0, 0, 0, 0) and are discussed in detail in section 5.2 in the context of the study of double instability.
If we compute a larger number of intersections we observe a drastic change in the topology of the ribbon-like surface. The points deviate from the well-defined 4D structures that we described above and form a cloud of points around them in the 3D projections with mixed colors in their 4D representations. However, in the particular case we study here, for about 3580 intersections, i.e. for 1000 consequents more than those building the structures in Figs. 2 and 3, in the 3D projections that include the z andż coordinates we discern a Θ structure (i.e. a greek letter Θ rotated by 90 o ) surrounded by a ring. This can be seen in Fig. 4, where we depict the consequents of the orbit "u" in a (ẋ, z,ż) projection and color them according to their x values. This new structure in the 3D projection of the phase space is maintained for several thousands of intersections more. We have again a confinement of the consequents in a restricted volume of the phase space, which is larger than previously. After 3580 intersections the consequents extend within a region (∆x × ∆ẋ × ∆z × ∆ż) < (0.18 × 0.7 × 0.8 × 2.4). In the new situation we have on one hand the appearance of structures in certain projections which indicates order, and on the other hand mixing of colors, which indicates chaoticity.
The origin of the observed structure can be understood best in the (z,ż) projection (Fig. 5). We remind that we study the structure of phase space in the region U, which starts with the bifurcation of the family x1v2 at the U → S transition of x1 in the vertical 2:1 resonance region (Fig. 1). For smaller energies in the resonance region, at E j = −5.028, we have the S → U transition of x1 and the introduction in the system of the stable family x1v1 (and its symmetric family x1v1 ′ ). Thus for E j = −4.66, we have, apart from the two simple unstable p.o. x1v2 and x1v2 ′ , the two stable x1v1 and x1v1 ′ and the planar stable x1. The presence of all these simple periodic orbits determines the orbital behavior of "u" in the following way: In the central part of Fig. 5a we observe a red "ellipse with corners", which forms within the first 2580 intersections and which is the projection of the structure of Fig; 3 in the (z,ż) plane. The initial conditions of x1v2 and x1v2 ′ are located at the upper and lower corner of the red elliptical curve respectively. The blue ellipse inside the red one is the projection of a rotational torus (Vrahatis et al. 1997) belonging to a quasiperiodic orbit found by perturbing x1 in theż direction. During the first 2580 consequents the orbit "u" stays sticky to the outer 4D rotational tori of x1. The outermost rotational tori belonging to x1 are broader than the ribbon-like tori of the "u" orbit. In Fig. 6 we give two characteristic projections from the point of view (θ, φ) = (90 o , 42 o ) in spherical coordinates. In (a) we have the (x,ẋ, z) and in (b) the (ẋ, z,ż) projection. In both of them we overplot the 2580 consequents of "u" (red) and one of the outermost rotational tori around x1 (black).
By continuing the integration of the orbit "u" for more than 2580 intersections we see that the next 550 intersections form a Θ -like structure with two lobes on the sides of the red elliptical structure (Fig. 5a, indicated with a couple of white arrows). These consequents surround the rotational tori around the stable (a) 5000 consequents of the orbit "u" in the 2D (z,ż) projection. The first 2580 consequents form the red "ellipse with corners" in the central part of the figure. The location of x1v2 and x1v2 ′ are indicated with black symbols on it (up and down). The blue ellipse inside the red one is the projection of a rotational torus around x1. The rest of the consequents, after the first 2580, diffuse in phase space, but remain in the region around the rotational tori belonging to the stable p.o. of the x1v1 and x1v1 ′ families. The projections of the two tori around each of the two stable p.o. (indicated also with black symbols in the middle of the "empty" regions) appear in this case like green and magenta invariant curves. Arrows point to the lobes and ring structure (see text). (b) 15000 consequents of the orbit "u" in the 2D (z,ż) projection. They visit all available regions of the phase space. However the sticky regions formed during the integration time are clearly discernible.
p.o. x1v1 and x1v1 ′ , which are very thin and appear in this projection like invariant curves (green and magenta). The rotational tori of these two families are still in the phase space region close to the U p.o. x1v2, since they have been bifurcated at a nearby E j from x1 (Fig. 1). The initial conditions of x1v1 and x1v1 ′ are in the middle of the rotational tori, marked with black symbols.
Considering more than 3130 consequents we find that the next 1370 intersections stay trapped in a ring in the (z,ż) projection, which in its turn surrounds all structures we have found in the phase space up to now. We point to this ring with two black arrows in Fig. 5a (right above and left below). Between the two lobes surrounding the rotational tori around x1v1 and x1v1 ′ and the ring structure, we observe a chain of eight stability islands (colored black), which are the (z,ż) projections of the corresponding rotational tori of a stable 8-periodic orbit. Finally, after 4500 intersections, the consequents start diffusing into a larger volume of the phase space. Figure 5a depicts the space of section after 5000 intersections. After 15000 intersections the consequents have occupied all available phase space (Fig. 5b). However, the denser regions clearly show us where the orbits have been trapped for longer times.
We have to compare the structures described up to now with the 4D space of section we presented in Fig. 4. The projections in Fig. 5 are part of Fig. 4. We observe them as we rotate the (ẋ, z,ż) diagram in our computer screen. Since Fig. 4 consists of 3580 consequents, apart from the ribbon it includes also the consequents forming the Θ -like structure around the rotational tori of x1v1 and x1v1 ′ as well as part of the consequents that form the "ring" in Fig. 5. In Fig. 4 we emphasize on the distribution of the consequents in the 4th dimension. We observe that in Fig. 4 the color, representing the normalized x-coordinate values, is dominated by green shades. This means that although we do not have a smooth color variation as in Figs. 2,3 the consequents stay confined in a restricted volume of phase space. This happens indeed. The x values for the first 3130 intersections are in the range −0.008 < x < 0.08, while within the next 450 intersections this range is almost doubled. Due to this expansion of the range of x values, in Fig. 4 are introduced the orange-red and blue-violet consequents. The fact that in 3D projections we find structure and that the color variation is dominated by a particular shade, indicates a weakly chaotic behavior of the orbit "u". This information is included in the 4D representation of the phase space as it is given in Fig. 4. By combining Figs. 4 and 5 on a computer screen, we reveal not only the weakly chaotic character of the orbit, but also the regions where it sticks during the integration time. This shows the efficiency of the color-rotation method.
One can calculate for the orbit "u" the finite-time Lyapunov Characteristic Number, LCN (t), which is defined by where ξ(t 0 ) and ξ(t) are the distances between two points of two nearby orbits at times t = 0 and t respectively (see e.g. Skokos 2010). As we can see in Fig. 7 the changes we have found in the orbital behavior of "u" do correspond to characteristic changes of the slope of the LCN (t) curve. During the first 2580 intersections, indicated by the left arrow in Fig. 7a the curve decreases towards zero. Then during the next 450 intersections, until the time pointed by the right arrow in Fig. 7a when the consequents form the Θ -like structure with the lobes around the x1v1 and x1v1 ′ tori, it continues to decrease but the curve has a much smaller slope. After that it starts increasing and its evolution after this point until t = 3 × 10 4 can be followed in Fig. 7b (note the different scales in the vertical axes). LCN (t) increases, then levels off, and has a decreasing part until the left arrow in Fig. 7b. which corresponds to about 4180 intersections. At this point we have an expansion of the consequents in phase space, which can be understood in the (x,ẋ) projection (not shown in a figure). Beyond the left arrow in Fig. 7b the LCN (t) curve increases. Then, between Fig. 7b and 7c the curve decreases and in Fig. 7c it increases again. Finally, for longer times of integration, when the orbit explores all the available phase space, the curve levels off at a value LCN (t) ≈ 2.74 × 10 −2 (Fig. 7c) The right arrow in Fig. 7b corresponds to 4500 intersections and the region to the right of this arrow corresponds to the expansion of the consequents from the ring to the whole available volume of phase space as we observe in Fig. 5b. Although LCN (t) gives globally the variation of the volumes of the phase-space filled by the consequents, by means of the color-rotation representation we can see in detail the specific regions where the consequents are trapped and the coordinate along which they diffuse to occupy larger volumes of the phase space. This information is essential in Galactic Dynamics for building self-consistent models using the orbital theory.
Summarizing the investigation of the orbital behavior of "u" in the neighborhood of a simple unstable periodic orbit we underline that it is characterized by two sticky periods, one associated with the rotational tori of x1 and a second sticky period associated with the rotational tori of x1v1 and x1v1 ′ . During the latter sticky period we observe the Θ structure and the ring in the (z,ż) projection. In the following section we examine how general is this behavior.
Another class of orbits in the neighborhood of x1v2 is characterized by immediate stickiness around the x1v1 and x1v1 ′ rotational tori, without building initially a "ribbon". We find them by keeping the same energy constant (E j = −4.66) and computing orbits by adding the perturbations given in Table 1 to the initial conditions of x1v2. In all these cases we do not observe the well defined ribbon-like surface that we have seen in Figs. 2 and 3. The orbits, starting close to the p.o. x1v2, depart directly from its immediate neighborhood and stick around the rotational tori of x1v1 and x1v1 ′ , before they expand and explore all Table 1. The range of deviations from the initial conditions of the periodic orbit x1v2 for E j = −4.66 that give orbits whose consequents stay sticky to x1v1 and x1v1 ′ and then diffuse occupying all the available phase space.
available phase space. They form Θ -shaped structures surrounded by rings as in Fig. 4 without the central elliptical structure we observe in Fig. 5. In Table 1 we give the ranges of the perturbations of the initial conditions for which we have observed this orbital behavior. For deviations from the p.o. in theż direction in the interval 4.4 × 10 −1 < ∆ż < 4.7 × 10 −1 we reach the eight tori of the stability islands (Fig. 5). By perturbing an initial condition we keep all others equal to the corresponding initial condition of the p.o.
For perturbations larger than those in Table 1 in the x,ẋ orż-direction we do not find any stickiness but the consequents form clouds around the initial conditions of the p.o. x1v2. These clouds are not confined in a volume smaller than the total available phase space. They occupy always larger volumes as the integration time increases. We have to take into account that for larger ∆z perturbations than the listed in Table 1 we first fall on initial conditions of the rotational tori of x1v1 and x1v1 ′ and then for even larger perturbations we find again the chaotic sea.
The second way we followed to explore the structure of phase space around x1v2, is by varying E j , keeping the same perturbation in its initial conditions. We did so by adding a ∆x = 10 −4 perturbation in its x 0 initial conditions, for E j 's in the region U. For −4.7 < E j < −4.4 the orbital dynamics of the orbits are as in the case for E j = −4.66 that we described in section 4.1. However, for E j > −4.4 a cloud of points surrounds the initial conditions of x1v2 in the 4D phase space. This means that in all 3D projections we observe clouds surrounding x1v2 from the beginning of the calculation and also a mixing of colors.

Asymptotic Curves
Let us now try to understand the orbital behavior of the orbits we have studied in the previous section by investigating the structure of the asymptotic curves of x1v2 in the 4D spaces of section. The asymptotic curves are computed by mapping a line segment in the space of section on an associated dilating eigenvector, starting very close to a simple unstable periodic orbit x1v2 (see e.g. Magnenat 1981). The example we present is again for E j = −4.66. We have already studied the structure of the phase space in the neighborhood of this simple unstable periodic orbit in sections 4.1 and 4.2. The initial conditions of the p.o. have been given in section 4.1.
We consider the unstable eigenvector (x 1 ,ẋ 1 , z 1 ,ż 1 ) = (−0.0011989, −0.00109617, −0.0459267, 0.998943) and the corresponding eigenvalue λ = 1.63536. We compute the asymptotic curves of the unstable manifold and we study the morphology of its asymptotic curves. The initial segment is composed of the initial conditions (x 2 ,ẋ 2 , z 2 ,ż 2 ) of the form: First we compute the asymptotic curves for c < 0 and we call this part of the unstable manifold, part A. The initial segment for part A is composed of 10 4 points (for −10 −3 ≤ c ≤ −10 −7 with a step 10 −7 ). For each point, 50 consequents were calculated. The basic criterion for choosing the values of the perturbations and the number of the consequents is to stay on the asymptotic curves with sufficient accuracy. The monodromy matrix conserves its symplectic character (determinant = 1) with an accuracy at least of the order of 10 −6 .
In Fig. 8 we observe the first 50 consequents of part A of the unstable manifold in the (x,ẋ, z) subspace of the 4D space of section. The first 2 consequents are very close to the periodic orbit. Then, we discern the 3rd consequent (indicated with an arrow). The next 13 consequents follow the direction of the black arrows, starting to the right and upwards, and they form a loop. After that, we return back to the neighborhood of the periodic orbit. Then, the 16th consequent as well as the following 7 remain again very close to the periodic orbit. From the 24th consequent and for the next 14 consequents we observe the formation of another loop, below the first one, following the directions that are indicated by the arrows. The final 12 consequents are again very close to the periodic orbit. The evolution in time of an initial segment of initial conditions leads to the formation of a ribbon-like double loop structure, similar to the one in Fig.  2. The same holds for the (ẋ, z,ż) 3D projection of the manifold, which has a similar structure with the corresponding projection of the orbit "u" (Fig. 3). As in Figs. 2,3 color on the asymptotic curves varies smoothly, so the 4D morphology is as for the orbit "u".
For c > 0 in the initial segment of the initial conditions the asymptotic curves evolve in the way we describe below and form the second part of the unstable manifold. We call this second part, part B. The initial segment of initial conditions for part B is composed by 2 × 10 4 points (for c = 10 −6 ≤ c ≤ 2 × 10 −2 with a step 10 −6 ). For every point we computed this time 20 consequents. In Fig. 9 we depict the first 60 points of the initial segment in the (x,ẋ, z) projection and color them according to theirż values. The figure helps us understand that the asymptotic curves of part B of the unstable manifold are 4D curves, which wind around the point corresponding to the x1v2 p.o. in the space of section. Note that we have a smooth color variation along the curve from light blue to red, following the direction of the successive numbered arrows.
Evolving the orbits for a larger number of points (consequents) from the initial segment (400 instead of 60 points) we observe that the asymptotic curves of part B oscillate rather far in the (x,ẋ, z) projection (Fig. 10). For example in the z-direction these oscillations extend between −0.4 < z < 0.4. Considering the evolution of all 2 × 10 4 consequents of the initial segment the asymptotic curves fill roughly the 4D space (x 1 , x 2 ) × (ẋ 1 ,ẋ 2 ) × (z 1 , z 2 ) × (ż 1 ,ż 2 ) = (0.025, 0.05) × (−0.08, 0.08) × (−0.4, 0.4) × (−1, 1) (Fig. 10). The figure is complicated and one can discern apart from the formation of loops at the upper left corner in Fig. 10, a ring/disk structure that has at its center the initial conditions of the p.o.  As in the case of the orbits we have discussed in subsections 4.1 and 4.2, we also study the asymptotic curves by combining several 3D projections in which we color their points according to the fourth coordinate. Using the (ẋ, z,ż) 3D projection and coloring the points according to their x values the asymptotic curves are depicted in the following way.
In Fig. 11 we give the first 400 consequents totally (20 points on the initial segment evolved for 20 intersections). They form an asymptotic curve that oscillates in the 3D (ẋ, z,ż) projection. Along the asymptotic curve the color variation is smooth (from light blue to red to dark blue and then back to light  blue), which means that the asymptotic curve is a 4D curve. Evolving further the asymptotic curve in time we observe that the oscillations become more complicated but the smooth color variation continues along the curve. This clearly shows that the oscillations, which have been calculated by Magnenat (Magnenat 1982) in 2D projections, occur in the 4D space.
In order to study the large scale evolution of part B of the unstable manifold we evolve next 400 points of the initial segment. The result is depicted in Fig. 12. After one oscillation on the left side of Fig. 12a we have several oscillations in the right part of the figure. The oscillations become longer than those of the previous case (Fig. 11) and they extend to larger absolute values ofż in the 3D subspace (ẋ, z,ż) (Fig.  12a). By evolving 20000 points of the initial segment, the asymptotic curves form a Θ -shaped surface in the (ẋ, z,ż) projection (Fig. 12b) consisting of many parallel filaments. The filamentary character of this structure is due to the parallel oscillations of the asymptotic curves in the (ẋ, z,ż) subspace.
There is an obvious similarity of Fig. 12 with Figs. 4 and 5. The projection of the asymptotic curves in five of the six 2D projections show clouds of points that seem to fill densely the available phase space. However, in the (z,ż) projection we get a clear insight of how the curves evolve in time. In this projection (Fig. 13) we see that the asymptotic curves of part B of the unstable manifold depart from the p.o. in the unstable direction (arrow 1) and then they follow the arcs as indicated by the arrows 2,3,4,5 and 6. The asymptotic curves oscillate along a direction vertical to the stable eigendirection.
The morphology we encounter for the case of the stable manifold is similar to that of the asymptotic curves of the unstable invariant manifold. Below we discuss the morphology of the asymptotic curves for other values of the Jacobi constant in the region U of simple unstable x1v2 periodic orbits.

Energy variation and asymptotic curves
We describe now the invariant manifolds in the neighborhood of a simple unstable periodic orbit for values of E j −4.7 < E j < −4.62. They have again two parts.
As the energy increases, for values larger than E j = −4.66, the asymptotic curves of part A form a ribbon, double loop, surface in the 3D subspace (x,ẋ, z). As an example we give the case with E j = −4.63 (Fig. 14). Although the lobes of the double loop are distorted, they form an acute angle in their (x, z) projection like the corresponding lobes of the orbit "u" (lower middle panel in Fig. 2). We observe also that the consequents wind around the periodic orbit (as the consequents in Fig. 9) when they approach the periodic orbit. The manifold has a similar structure like the one we have encountered in the central region of "u" in the same 3D projection. Also, in these higher energies, the asymptotic curves of part A form again a circular ribbon-like surface in the 3D subspace (ẋ, z,ż) similar to the orbit in Fig. 3. The increase of energy causes oscillations of the asymptotic curves of part A, in the 3D subspace (ẋ, z,ż).

Stickiness
From the investigation of the spaces of section and the numerical calculation of the asymptotic curves we presented in the two previous subsections it becomes evident that the phenomenon of stickiness not only is present but characterizes the dynamical behavior of many orbits that start in the neighborhood of the simple unstable orbit x1v2. As we know, stickiness in chaos refers to chaotic orbits that stay in a particular region of the phase space, before escaping to large distances. In 2D autonomous Hamiltonian systems stickiness appears near the borders of stability islands, as well as near the asymptotic curves of the unstable manifolds for a long time (Contopoulos & Harsoula 2008, Contopoulos & Harsoula 2010). In our paper we explore this phenomenon in the 4D space of section of our 3D system and we try to explain it  Fig. 12b. The arrows show the direction followed by the consequents. The stable and unstable eigendirections are depicted with magenta and a green lines respectively. The periodic orbit x1v2 is indicated with a black "+" symbol. in combination with the structure of the associated manifolds. It is the first time that stickiness in chaos is studied in detail in the neighborhood of simple unstable periodic orbits of a 3D autonomous Hamiltonian system. In order to study the phenomenon of stickiness in our system we consider first the representative example of a periodic orbit in the region U of the x1v2 family for E j = −4.66.
If we add a negative perturbation in theż direction (as we mentioned in section 4.2), or add a small positive perturbation in theẋ direction (for 10 −6 ≤ ∆ẋ < 2 × 10 −4 ) to the initial conditions of the periodic orbit x1v2, the orbits we find have a dynamical behavior similar to the one of the orbit "u" (section 3.1). By comparing the orbital behavior of all orbits having a morphology like "u", with the form of the unstable manifolds we conclude that the initial conditions of "u" are very close to the unstable eigendirection of the manifold. Both the consequents of "u" and the asymptotic curves wind around the rotational tori associated with x1. The consequents of "u" follow the asymptotic curves and stick in the region dictated by the manifolds. For example in Figs. 15a and 15b the first 2580 consequents of the orbit "u" (depicted by blue color) stick close to the asymptotic curves of part A of the unstable manifold (drawn with red color). This number of consequents correspond to a time of 124 periods of rotation of our system, a very large time interval in galactic dynamics (of the order of several Hubble times).
For a larger number of consequents, the points depart from part A of the unstable manifold. The relation between the consequents of "u" and the unstable manifold is best understood again in the 3D projection (ẋ, z,ż). As we can see in this projection in Fig. 16, the first 1000 consequents, after the departure of the consequents from part A of the unstable manifold, remain close to the asymptotic curves of part B of the unstable manifold. This is the mechanism that creates in phase space the lobes of the Θ and the ring structure we observe as denser regions in the (z,ż) projection in Fig. 5. These denser regions are created due to the phenomenon of stickiness around the rotational tori of x1v1 and x1v1 ′ . The lobes of the Θ structure, indicated with white arrows in Fig. 5, are formed by consequents that stay almost on the (red) asymptotic curves of Fig. 16. However, even after the departure of the consequents from the lobes they stay close to the x1v1 and x1v1 ′ rotational tori, because they are trapped in the immediate neighborhood of the eight islands of stability we observe as small black line segments in Fig. 5. In this projection the figure resembles what is found in 2D systems close to the last KAM curve (Contopoulos 2002). However, we should have in mind that in the 4D space of section of our 3D Hamiltonian system consequents can be found on both sides of a torus, moving through the additional dimensions.
In section 4.2 we said that for the values of perturbation listed in Table 1, the orbits go directly around the rotational tori of x1v1 and x1v1 ′ . All these orbits have also a sticky behavior at the borders of the Fig. 16. The first 1000 consequents of the orbit "u", given with green "+" symbols after their departure from part A of the unstable manifold, as observed in the 3D projection (ẋ, z,ż). These consequents stick close to the asymptotic curves of part B of the unstable manifold (drawn with red color). Our point of view in spherical coordinates is (θ, φ) = (78 o , 79 o ). x1v1 and x1v1 ′ 4D stability islands, which we can always see best in the 3D projection (ẋ, z,ż) and the 2D projection (z,ż). In all other projections the sticky regions are discernible as denser regions of the phase space embedded in clouds. The (ẋ, z,ż) projection has the advantage that in a certain range of viewing angles only a few consequents from the cloud are projected in the regions, which look "empty" (actually they are occupied by the rotational tori around the stable orbits).
Our study of the dynamical behavior in the neighborhood of the simple unstable orbits, shows that the characteristic double loop structures that have been already found by Magnenat (1982) and Patsis and Zachilas (1994) etc., as well as the Θ -shaped structures and rings we observe in Fig. 5, are associated with stickiness at the borders of rotational tori of stable families that exist in the system. The consequents follow the manifolds that wind around the tori. Thus it would be more accurate to associate these structure of the phase space mainly with the U → S transitions, as they are not expected to exist away from such a transition point or after a DU → U transition.

Double Instability
Double instability in rotating galactic potentials usually appears for large energies. In the x1 orbital tree  DU appears in higher order bifurcations, when a stable 3D bifurcation of x1 becomes first simple unstable and then double unstable, or when a simple unstable bifurcation of x1 becomes double unstable. Frequently we find DU p.o. at energies for which nearby orbits escape. Thus the selection of the appropriate case for investigating the phase space structure is crucial. Here we consider first the dynamics in the neighborhood of orbits of the x1v2 family in the region DU (Fig. 1). However, in order to examine the dependence of the dynamical behavior on the energy of the orbit, we study also the phase space close double unstable z-axis orbits. These are p.o. along the rotation axis z in rotating triaxial potentials (Heisler et al. 1982).
Then we explore the structure of the phase space in the neighborhood of the double unstable periodic orbits of x1v2 in the region DU from −3.41 < E j < −2.3, for the same fixed value of the perturbation, ∆x = 10 −4 . We find again 4D clouds of points as in the case of the double unstable periodic orbit at E j = −3.38.
We study the morphology of the asymptotic curves of the unstable manifold in the neighborhood of the double unstable periodic orbits in the region DU and we present this morphology in the neighborhood of the double unstable periodic orbit for E j = −3.38 (Fig. 17). The computation of the asymptotic curves is similar to that of the case of the simple unstable periodic orbits (section 4.3). The only difference is that we have now two unstable eigenvectors. These are ( x 1 ,ẋ 1 , z 1 ,ż 1 ) = (0.0216618, 0.163481, −0.737238, 0.655198) and (x 2 ,ẋ 2 , z 2 ,ż 2 ) = (−0.831661, 0.523481, 0.0173103, 0.184412). The corresponding eigenvalues are λ 1 = 69.08801 and λ 2 = 1.74981. The initial conditions (x 3 ,ẋ 3 , z 3 ,ż 3 ), that we use for the computation of the asymptotic curves, are given by: We computed the asymptotic curves, in two cases which we named "I" and "II" respectively. In the first case, we kept c 2 constant and equal to zero and we used c 1 as a parameter. In the second case, we kept c 1 = 0 and we varied c 2 . In case "I" we have taken values c 2 = 0 and 10 −7 ≤ c 1 ≤ 2 × 10 −4 with a step 10 −7 (2000 initial conditions). In the second case "II" we took values c 1 = 0 and we varied c 2 in the interval −1.35 × 10 −1 ≤ c 2 ≤ −10 −5 with a step 10 −5 (13500 initial conditions). In both cases "I" and "II" we computed only 3 consequents for every initial condition 1 . Here we present the morphology of the asymptotic curves in the 3D (x,ẋ, z) projection. We found in our numerical experiments, that the morphology of the asymptotic curves in all 3D projections is similar. As it was expected the consequents along the asymptotic curves in case "I" (with the large eigenvalue λ 1 ) move away faster from the periodic orbit than the consequents along the asymptotic curves in case "II" (with the small eigenvalue λ 2 ).
The consequents of the unstable invariant manifold in case "I" depart from the periodic orbit by following the black arrows of Fig. 17. They form first a loop in the 3D projections of the 4D space of section below the periodic orbit. Then the consequents approach again the periodic orbit but with a different color (yellow) than the one of the initial part which is violet. Afterwards the asymptotic curve performs a larger loop before leaving the box of the frame. In addition we observe a smooth color distribution from violet (close to the p.o.) to blue, green, yellow and red and then follows a color oscillation from red to green and back again. These oscillations correspond to oscillations ofż in the upper part of the color spectrum (right part of Fig. 17). This smooth color variation along the curve indicates that these asymptotic curves are 4-dimensional curves.
In both cases "I" and "II" the asymptotic curves deviate soon from the eigendirections. In Fig. 18 we compare the asymptotic curves "I" (blue color) and "II" (red color) that start along the unstable eigendirections that correspond to the large (λ 1 ) and small (λ 2 ) eigenvalues respectively. The morphologies of these curves are different. The orbits we start integrating close to these curves remain sticky to them only for a short time (0.4 periods of rotation of our system) before moving away from them. Soon they fill in a chaotic way a large region of the phase space. We have found similar morphologies for all cases of stable and unstable manifolds we have examined in the DU region. This is a dynamical behavior typical for DU periodic orbits of the x1 tree. However, it is not necessarily typical for the dynamics in the neighborhood of DU periodic orbits in general. We investigated the orbital behavior close to double unstable p.o. in several cases in our system. Below we describe the dynamics close to z-axis DU p.o.'s, which differs in many aspects from what we have found in the neighborhood of the DU x1v2 p.o.'s we presented up to now.  x .
z Fig. 18. The 3D subspace (x,ẋ, z) of the asymptotic curves of the unstable invariant manifold for case "I" (blue color) and case "II" (red color), in the neighborhood of the double unstable periodic orbit for E j = −3.38 in the 4D space of section. The two unstable eigendirections that correspond cases "I" and "II" are depicted with a magenta and a green line respectively. The location of the periodic orbit is given with a black square symbol. The point of view in spherical coordinates is (θ, φ) = (63 o , 168 o ).

z-axis
As already mentioned, the orbits of the z-axis family lie entirely on the rotational, z axis, of the system. Thus, we consider now the z = 0 surface as surface of section withż > 0. Earlier studies have shown that this family becomes double unstable in slowly rotating triaxal systems (Martinet & de Zeeuw 1988, Patsis b Fig. 19. (a) The 2D (x,ẋ) projection of the orbit we study in the neighborhood of the DU z-axis orbit at E j = −4.9. We consider here 4200 consequents. The first 250 are plotted with blue color. The next 3950 form the red "8"-shaped structure.
Then they start to diffuse in phase space. Black symbols indicate the location of the periodic orbits in this projection. (b) The first 4200 consequents are depicted now with heavy red points, while the next 26000 (light red dots) diffuse and occupy all available phase space.
& Zachilas 1990). We investigated the stability of this family in the potential (2) rotating very slowly with Ω b = 5 × 10 −6 . For this value of Ω b the z-axis family has large regions of double instability (DU). In this case, as we vary the energy, the z-axis family for large negative values of E j is stable and for E j = −6.7 it becomes simple unstable. At this transition bifurcates the stable family of p.o. of stable anomalous orbits, abbreviated to sao (Heisler et al 1982), and its symmetric family with respect the yz-plane. Then the z-axis family, for E j = −5.836, becomes double unstable and it remains DU for larger values of E j . In the (x,ẋ, y,ẏ) 4D space of section the initial conditions of the z-axis orbits are (0, 0, 0, 0). For E j = −4.9 we present the case of the orbit in the neighborhood of the z-axis p.o., which deviates from it in theẋ direction by ∆ẋ = 0.62. We call this orbit "duz1".
In the 2D (x,ẋ) projection we observe the first 250 consequents being distributed around the p.o. (blue points in Fig. 19a). These points do not show clearly any particular structure. The next 3950 consequents however, form an "8"-figure. We plot with red color 950 of them in Fig. 19a. In the same figure the periodic orbits existing in this energy are given with black symbols. In the center of the "8"-figure structure the × symbol indicates the DU z-axis orbit, while in the central regions of the two lobes of this structure the black dots (on the left and the right) indicate the location of sao and its symmetric p.o., which are surrounded by projections of rotational tori around them. If we continue to integrate the orbit, we find that after more than 4200 intersections the consequents diffuse and occupy all available phase space. Then the (x,ẋ) projection is as we see in Fig. 19b.  Table  2. The range of perturbations that give orbits whose consequents are distributed as in Fig. 19.
On one hand we have the formation of a structure discernible in projections of the phase space. On the other hand, the consequents that form the "8"-figure structure of Fig. 19a are not on a smooth 4D surface. The application of the color-rotation method shows that we can hardly observe a smooth color variation along any direction. These are indications of a chaotic behavior.
The finite time LCN index points also to the chaotic character of the orbit. For most of the time when the consequents in the (x,ẋ) (or the (x,ẋ, y)) projection stay on the "8"-shaped figure, the index varies around 5 × 10 −2 and finally decreases until the time pointed by an arrow in Fig. 20a. This time corresponds to 4200 consequents. In Fig. 20b we see how LCN (t) evolves as the consequents start filling all available phase space. It increases continuously up to t = 28500 and beyond.
The range of perturbations for which we find the "8"-figure is registered in Table 2. Outside the range given for the x coordinate we find orbits sticky to rotational tori around the sao and its symmetric orbit. For absolutely larger perturbations in theẋ, y orẏ-direction we encounter three types of orbits, namely (a) orbits represented by clouds of points, (b) orbits sticky to tori around the stable periodic orbit x1v1 and (c) orbits sticky to chains of rotational tori belonging to p.o. of higher multiplicity. This orbital behavior is related with the unstable manifold of the z-axis orbit as follows.
The computation of the manifolds is similar to that of the case of the double unstable periodic orbits of x1v2. The unstable eigenvectors are (x 1 ,ẋ 1 , y 1 ,ẏ 1 ) = (−0.20427, 0.000957614, 0.978912, −0.00197504) and (x 2 ,ẋ 2 , y 2 ,ẏ 2 ) = (−0.000880093, 0.436429, 0.00422376, −0.899728). The corresponding eigenvalues are λ 1 = 5.93858 and λ 2 = 5.90405. The initial conditions (x 3 ,ẋ 3 , y 3 ,ẏ 3 ), which we use for the computation of the asymptotic curves, are given by: We calculate the manifold by using for c 1 the values −10 −3 ≤ c 1 ≤ 10 −3 with a step 10 −5 and for every value of c 1 we consider for c 2 the values −10 −3 ≤ c 2 ≤ 10 −3 again with a step 10 −5 . For every initial condition we compute 10 consequents. Then we have totally 4 × 10 5 consequents. This way we computed the unstable asymptotic surface we present in Fig. 21. We note that now we have an asymptotic surface and not an asymptotic curve. This is expected since we have two eigenvalues outside the unit circle (see e.g. Arnold 1988, p. 287). The asymptotic surface is easier to be visualized than in the case of the orbit in the neighborhood of the double unstable x1v2 orbit, because now the two associated eigenvalues are very close. We observe that the unstable asymptotic surface we calculated has a morphology of an "8"-figure in the 3D projection (x,ẋ, y) of the surface of section. For the 4th dimension we color the points according to the normalizedż values. We have a complicated but smooth color variation on the surface of the "8"-shaped structure (Fig. 21a). In Fig. 21b we give the (x,ẋ) projection of the same figure that shows clearly the central oscillations of the surface close to the p.o. located at the center of this "8"-figure structure.
A similar morphology is found for the stable asymptotic surface. The structure we present in Fig. 21 is part of the unstable manifold. If we consider a larger number of consequents for −1 ≦ c1, c2 ≦ 1 and a step 10 −3 in (5), the asymptotic surface expands considerably also in the y-direction and its (x,ẋ) projection is as the red points show in Fig. 22a. We realize that all consequents of "duz1" (black color) stick on the unstable asymptotic surface.
We have to note also that inside the lobes of the "8"-shaped manifold structure, besides the sao rotational tori (green elliptical curve in Fig. 22b) we find also sticky orbits, initially forming toroidal objects in this region and then diffusing to a larger volume, remaining however sticky to the surface of the manifold. The consequents of two such orbits are given with blue and magenta colors in the right lobe of the manifold in Fig. 22b. The existence of the asymptotic surface decelerates the diffusion of these orbits to larger phase space regions for long times.

Diffusion in phase space
The study of the orbits in the neighborhood of simple and double unstable periodic orbits gives us an information of practical interest namely their rate of diffusion and the possible trapping of their consequents in particular regions of the available phase space. This would allow us to separate the orbits that support physical structures from those which do not. For this reason we measure the diffusion of the orbits we study and we compare the results among themselves.
From the orbits close to x1v2 p.o. we have consider the following orbital behaviors: • The first type of orbits we found close to x1v2 periodic orbits is represented by double loop ribbons or The (x,ẋ) 2D projection of the unstable asymptotic surface (red color) and the consequents of orbit "duz1" that remains for long times sticky to it. (b) The same projection of the unstable asymptotic surface, together with two orbits (blue and magenta consequents) that we find by perturbing orbits of the sao rotational tori (projected to a small green elliptical curve near the center of the right part of the manifold).
by closed ribbons in the 3D projections of the 4D space of section. An example is the orbit "u" presented in section 4.1. (Figs. 2,3). The orbits of this orbital type continue beyond the initial ribbons, by forming additional dense regions in phase space after long integration times. These dense regions are discernible in some projections as "8"-like structures or rings (Fig. 5). • The second type of orbital behavior is represented by clouds of points in the 4D surface of section. As an example, we consider the orbit with ∆x = 10 −4 from the initial conditions of the simple unstable x1v2 for E j = −4 and we call it "u1". • The third type of orbital behavior is represented again by clouds of points in the 4D surface of section, but this time in the neighborhood of a double unstable periodic orbit. As a typical example for the discussion to follow, we consider the orbit deviating from x1v2 by ∆x = 10 −4 at E j = −3.38. We call this orbit "du1".
In addition we include an orbit in the neighborhood of a z-axis DU p.o., i.e.
• The orbit that builds an "8"-shaped figure before diffusing and occupying all available volume of the phase space (Fig. 19), which we call "duz1".
In order to estimate the diffusion velocity for the above types of orbital behavior we define a diffusion time and a diffusion velocity as follows: In general the volume of phase space occupied by the consequents of a chaotic orbit increases with time until they occupy the maximum available volume restricted by the surface of the zero velocity. The time that the orbits need to occupy the maximum volume V max , we define as "diffusion time" T d . Then the "diffusion speed" is the ratio The quantities V max and T d for the orbits in the neighborhood of a periodic orbit can be easily computed in two steps: (1) As time increases we compute the mean distance from the periodic orbit s up to this time t of all consequents and then we compute a mean volume V (t). In other words we compute a "mean volume" of the hypersphere centered on the periodic orbit with radius the mean distance s, i.e the quantity: V (t) = 1 2 π 2 s 4 . (2) After time T d the consequents occupy the maximum volume V max . This is represented by the maximum peak in a V (t) diagram and this allows us to compute the diffusion speed.
In Fig. 23a we plot the volume occupied by the consequents of the orbit "u" (Figs. 2,3) versus time t for an interval t ≤ 28000, corresponding to 4200 intersections. We observe that for t 17300 in our system  Fig. 23a corresponding to 2580 intersections) V (t) remains close to 0.165 after some initial oscillations. This time interval corresponds to the phase that the consequents stick to the asymptotic curves of part A of the unstable manifold forming the double loop ribbon (section 4.1). Then we observe a second interval between the two arrows in Fig. 23a, where V (t) increases but the slope of the curve is smaller than the slope beyond the right arrow. The second arrow corresponds to time t = 21800 or 3130 intersections. The dynamical behavior between the two arrows corresponds to the time during which the consequents form the lobes of the Θ -structure and the "ring" in the (z,ż) projection (Fig. 5) sticking on the asymptotic curves of part B of the unstable manifold (section 4.4). Beyond the right arrow V (t) increases continuously with a characteristically steeper slope than before. The consequents leave all structures observed in the neighborhood of x1v2 (section 4.1). Figure 23b is a continuation of Fig. 23a and shows that the peak of V (t) corresponds to V max = 13.433 for T d = 510000. The speed of diffusion is u d = Vmax T d = 2.63 × 10 −5 . After the peak of V (t) in Fig. 23b the consequents arrive at distances from the periodic orbit that are smaller or equal to the maximum distance. This means that V (t) decreases a little and finally, it levels off to a constant value. This happens for all chaotic orbits.

units (indicated by the left arrow in
The evolution of V (t) and the quantity T d change in the case of the orbit "u1". In Fig. 24 we see that for the curve labeled "1", corresponding to the orbit "u1", we have V max = 117.397 and T d = 3180. This means that the speed of diffusion is now u d = 0.0369, i.e. it is 1403 times larger than the velocity of diffusion of orbit "u". The upper, red curve, labeled "2", describes the evolution of V (t) for the orbit "du1" in the neighborhood of the DU p.o. x1v2 at E j = −3.38. We observe that the curve has a maximum Fig. 24. The evolution of V (t) for the orbit "u1" in the neighborhood of the simple unstable periodic orbit x1v2 for E j = −4 (1) and for the orbit "du1" (2), when x1v2 became double unstable for E j = −3.38.
V max = 458 and T d = 3400. In this case the speed of diffusion is u d = 0.135 that is 3.65 times and 5133 times larger than the diffusion speeds of the orbits "u1" and "u" respectively. For "u1" and "du1" V (t) increases rather fast to the peak of the diagram, indicating that we do not have any important stickiness in these cases.
Finally for the case we have found stickiness in chaos close to a DU z-axis p.o., i.e. for the orbit "duz1", the evolution of the quantity V (t) is given in Fig. 25. Initially V (t) decreases to a value about 1.68 during the first 1600 consequents. This period ends at the time indicated by arrow "1". This happens because the majority of the points occupy the right lobe of the "8" structure, as in the (x,ẋ) 2D projection (Fig. 19). This implies that the majority of the points are included in a volume smaller than the volume occupied by the consequents forming the whole "8"-shaped structure. Thus, the mean volume decreases. Then the next consequents complete the "8" structure by populating both lobes with similar number of consequents, so V (t) increases until V (t) = 3.11969, which corresponds to 3200 intersections (arrow numbered "2"). Then we observe a small plateau until the point indicated by the arrow labeled "3" (V (t) = 3.02259). The time at this point corresponds to 4200 consequents. This is the number of consequents that stay on the eight-figure surface before they depart from it and occupy a larger volume of the phase space. For larger times V (t) increases. In Fig. 25b we have a peak at V max = 15.77 for T d = 1.88 × 10 5 . The diffusion speed is u d = 8.385 × 10 −5 .
It is evident by comparing Figs. 23, 24 and 25, that the diffusion of chaotic orbits from the neighborhood of simple or double unstable p.o. does not depend on the kind of instability of the periodic orbit. The phenomenon of stickiness that decelerates the diffusion appears both to orbits starting close to simple as well as to double unstable periodic orbits. When this happens we observe plateaus in the V (t) variation.

Conclusions
In this paper we studied the structure of the phase space close to simple unstable and double unstable periodic orbits in a triaxial 3D system representing a rotating disk galaxy. We have encountered different orbital behaviors and different diffusion rates from the neighborhood of the periodic orbit. The phenomenon of stickiness seems to be ubiquitous and does not depend on the kind of instability of the periodic orbit in the center of the region we study. Our results refer to orbits in the neighborhood of the families x1v2  and z-axis (Martinet & de Zeeuw 1988) and can be summarized as follows: • The orbits starting in the neighborhood of a simple unstable periodic orbit form (a) double loop ribbons in the 4D space, (b) structures discernible in some projections after long integration times, (c) clouds of  (3)) winds around deformed rotational tori belonging to orbits deviating from the x1 orbits in the ±z or ±ż-direction. Part B (calculated for c > 0 in (3)) winds around the rotational tori of the symmetric x1v1 and x1v1 ′ p.o. • In both cases chaotic orbits are guided by the manifolds close to the regions occupied by the rotational tori associated with the stable periodic orbits existing in the region. These chaotic orbits are sticky and play a very important role, especially in Galactic Dynamics, since the stickiness time is in many cases of the order of a Hubble time or even larger. • Double loops and Θ -shaped structures surrounded by rings that are found in several 2D and 3D projections of the 4D spaces of section in 3D Hamiltonian systems are formed by orbits close to the asymptotic curves of simple unstable periodic orbits that surround the rotational tori of two nearby stable p.o. as in Fig. 5. Such configurations appear if the parent family (in our case the planar x1 family) has a S → U transition and a nearby U → S transition. At the first transition bifurcate two symmetric stable orbits (like x1v1 and x1v1 ′ ) and at the second transition we have the bifurcation of two simple unstable families (in our case the families x1v2 and x1v2 ′ ). Near the U → S transition of the parent family, the already bifurcated stable orbits are not far from the unstable orbit under study, thus the asymptotic curves of the latter families can surround their rotational tori. Away from the transition points of the varying parameter (in our case the Jacobi constant) the shape of the asymptotic curves become complicated and the consequents of the nearby orbits form clouds in the 4D spaces of section.
• Double instability in 3D rotating galactic disks appears for larger energies than simple instability in the families of the x1 tree. The usual case is after a U → DU transition and in general there are no rotational tori associated with stable families in the neighborhood of the DU p.o. Typically the consequents of the orbits form clouds of points in the 4D spaces of section. • In the case of a double unstable periodic orbit there are two different pairs of eigendirections, starting at the periodic orbit. The curves of each pair start in opposite directions. Between the eigendirections an "eigensurface" is formed. However the shape of this surface is complicated and it is difficult to be visualized when the associated eigenvalues are not close to each other. • In the case of z-axis DU periodic orbits in slow pattern speed models, at low energies, we find that the eigenvalues corresponding to the unstable eigenvectors are close to each other. Because of that the speeds of deviation along the two eigendirections are similar. This allows us to visualize the unstable eigensurface corresponding to the DU periodic orbit. We find sticky chaotic orbits whose consequents stay close to the unstable manifold for very long times (of the order of a Hubble time or more). We find these sticky orbits either by perturbing the initial conditions of the z-axis DU p.o., or by perturbing orbits on rotational tori associated with the stable anomalous orbit (sao) and its symmetric family. • Among the orbits we studied we found those close to the double unstable orbits of the x1v2 family having the largest diffusion speed. The sticky chaotic orbits close to the bifurcating point of the simple unstable x1v2 orbit and close to the double unstable z-axis orbit we have examined, have comparable diffusion speeds. They are much slower than the diffusion speeds of the orbits in the neighborhood of x1v2 simple unstable periodic orbits away from its bifurcation point, or of the double unstable orbits of the same family that have eigenvalues of very different size.
These results will be used in a forthcoming paper to assess the role of chaotic orbits in shaping the vertical profiles of disk galaxies.