Instabilities and stickiness in a 3D rotating galactic potential
Abstract
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 zaxis (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 zaxis 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. The sticky chaotic orbits close to the bifurcation point of the simple unstable x1v2 orbit and close to the double unstable zaxis orbit we have examined, have comparable diffusion speeds. These speeds are much slower than of the orbits in the neighborhood of x1v2 simple unstable periodic orbits away from the bifurcating point, or of the double unstable orbits of the same family having very different eigenvalues along the corresponding unstable eigendirections.
1 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 (Katsanikas et al. 2011). 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 (Patsis et al. 2002). The x1v1 family is introduced at the transition from stability to simple instability () as the stability of the central, planar family x1 undergoes a 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 transition, another family of periodic orbits (hereafter p.o.), x1v2, is bifurcated from x1 as simple unstable (Skokos et al. 2002). 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 zaxis 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 peanutshaped structure, which is characteristic of many disk galaxies observed edgeon (Patsis et al 2002). 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 3dimensional 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 and we consider intersections with . Hereafter by “intersections” we will refer to intersections, which fulfill these properties, until we start discussing the dynamics close to zaxis. Then we will change space of section. Every consequent in our cross section is determined by 4 Cartesian coordinates, e.g. . 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 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 .

A well defined structure in all 3D projections and a smooth color variation along them is evidence for regular motion.

Mixing of colors indicates chaotic behavior of the consequents in 4D (Patsis & Zachilas 1994, Katsanikas & Patsis 2010)
1.1 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 . and has 4 eigenvalues , , , , where , , with and . Thus, the stability of the p.o. depends on the two stability indices and and on the quantity .
For and , all four eigenvalues are on the unit circle in the complex plane and the orbit is called stable (see e.g. Contopoulos 2002, Contopoulos & Magnenat 1985).
Then we have three kinds of instability:

simple instability : and or , (two eigenvalues on the unit circle and two on the real axis)

double instability ): and and (four eigenvalues on the real axis)

complex instability () (all four complex eigenvalues are off the unit circle in the complex plane):
(for details see Contopoulos 2002, p.286). Hereafter we will indicate the periodic orbits as , , and according to their stability. A generalization of the three kinds of instability in Hamiltonian systems of 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 zaxis 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.
2 The Hamiltonian System
The Hamiltonian of our system in Cartesian coordinates is:
(1) 
where is the potential:
(2) 
This potential in its axisymmetric form can be considered as an approximation of the potential of the Milky Way (Miyamoto & Nagai 1975).
In our units, the distance =1 corresponds to 1 kpc. The velocity unit corresponds to 209.64 . The Jacobi constant corresponds to 43950 and we will call it “energy”. For the parameters we have used the following values: . For the angular velocity we used the value in our units. The parameters determine the geometry of the disks, while 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 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 nonescape orbits close to double unstable periodic orbits (Section 5).
3 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 and for successive periodic orbits in a family by changing only one parameter of the system (in our case ), while keeping the other parameters constant. In Fig. 1 we observe in such a diagram how the family of periodic orbits x1v2 (Skokos et al. 2002) bifurcates from the 2D central family x1 on the plane at the transition that happens in the vertical 2:1 resonance region (Contopoulos 2002 p.390). This bifurcation occurs for . We observe that the stability parameters , of the family x1v2 decrease as increases. The value of is initially about at (left side of diagram). When becomes we have a transition of x1v2 from simple to double instability () (for ) and the x1v2 periodic orbits become double unstable. The interval of x1v2, where it is simple unstable () and the interval of x1v2 where it is double unstable (), will be called hereafter regions U and DU respectively (Fig. 1).

4 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 and . However, we find qualitatively similar results by using the and 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 by varying the values and the direction (, , or ) of the perturbation. Secondly, we consider a constant perturbation and we apply this perturbation to the initial conditions of several periodic orbits of the x1v2 family in the region U for different values. We start by keeping constant and investigate as a sample case the dynamical behavior very close to the point where we have the transition of the x1 family and the bifurcation of the family x1v2 as .
4.1 Orbital behavior close to the bifurcating point  A sample case
The energy is indicated with an arrow in Fig. 1. The periodic orbit has initial conditions . In order to study the orbital behavior in its neighborhood we added a small perturbation in the direction equal to . 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 subspace (Fig. 2). In the bottom of the figure we give also from left to right the 2D projections and 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 in the projection (middle panel in the lower part of Fig. 2). The 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 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 . 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 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 ribbonlike 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 very close to the initial condition of the periodic orbit .
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 differ only in the sign of their initial conditions. In the 3D 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 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 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 form again a ribbonlike surface, however this time with a single loop morphology. Color is determined by the 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 space, since they have different initial conditions. Their positions are indicated with arrows in Fig. 3. The regions close to both of them are colored red (same 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 welldefined ribbonlike surface that has no selfintersections in 4D.
A similar analysis can be done by means of the two other 3D projections and , which give 4D structures morphologically similar to and respectively. It is evident that the orbit “u” in the neighborhood of a x1v2 p.o. for 2580 consequents is confined in a particular region of the space of section with dimensions (. 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 zaxis orbits in 4D spaces of section defined by . The zaxis orbits lie entirely on the rotational, z axis, of the system. They have initial conditions 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 ribbonlike surface. The points deviate from the welldefined 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 and coordinates we discern a structure (i.e. a greek letter rotated by ) surrounded by a ring. This can be seen in Fig. 4, where we depict the consequents of the orbit “u” in a projection and color them according to their 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 (. 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 () 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 transition of x1 in the vertical 2:1 resonance region (Fig. 1). For smaller energies in the resonance region, at , we have the transition of x1 and the introduction in the system of the stable family x1v1 (and its symmetric family x1v1). Thus for , 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 () 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 ribbonlike tori of the “u” orbit. In Fig. 6 we give two characteristic projections from the point of view in spherical coordinates. In (a) we have the and in (b) the 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 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 p.o. x1v2, since they have been bifurcated at a nearby 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 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 projections of the corresponding rotational tori of a stable 8periodic 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 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 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 values for the first 3130 intersections are in the range , while within the next 450 intersections this range is almost doubled. Due to this expansion of the range of values, in Fig. 4 are introduced the orangered and blueviolet 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 colorrotation method.
One can calculate for the orbit “u” the finitetime Lyapunov Characteristic Number, , which is defined by
where and are the distances between two points of two nearby orbits at times and 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 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 can be followed in Fig. 7b (note the different scales in the vertical axes). 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 projection (not shown in a figure). Beyond the left arrow in Fig. 7b the 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 (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 gives globally the variation of the volumes of the phasespace filled by the consequents, by means of the colorrotation 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 selfconsistent 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 projection. In the following section we examine how general is this behavior.
4.2 The general orbital behavior in region “U”.
A dynamical behavior similar to that of the orbit “u” is observed by perturbing the initial condition of x1v2 up to (i.e. an absolutely largest ). We find a similar behavior if we perturb by at most . Initial conditions that have the same , and with x1v2 and that deviates by , give rotational tori around x1 (see Katsanikas & Patsis 2010).
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 () 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 ribbonlike 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 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 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 , 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 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 , keeping the same perturbation in its initial conditions. We did so by adding a perturbation in its initial conditions, for ’s in the region U. For the orbital dynamics of the orbits are as in the case for that we described in section 4.1. However, for 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.
4.3 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 . 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
and the corresponding eigenvalue . 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 of the form:
(3) 
First we compute the asymptotic curves for and we call this part of the unstable manifold, part A. The initial segment for part A is composed of points (for with a step ). 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 ) with an accuracy at least of the order of .
In Fig. 8 we observe the first 50 consequents of part A of the

unstable manifold in the 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 ribbonlike double loop structure, similar to the one in Fig. 2. The same holds for the 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 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 points (for with a step ). For every point we computed this time 20 consequents. In Fig. 9 we depict the first 60 points of

the initial segment in the 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 projection (Fig. 10). For example in the direction these

oscillations extend between . Considering the evolution of all consequents of the initial segment the asymptotic curves fill roughly the 4D space (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 3D projection and coloring the points according to their 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 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 (Fig. 12a). By evolving 20000 points of the initial segment, the asymptotic curves form a shaped surface in the 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 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 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 . They have again two parts.
As the energy increases, for values larger than , the asymptotic curves of part A form a ribbon, double loop, surface in the 3D subspace . As an example we give the case with (Fig. 14). Although the lobes of the double loop are distorted, they form an acute angle in their 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 ribbonlike surface in the 3D subspace similar to the orbit in Fig. 3. The increase of energy causes oscillations of the asymptotic curves of part A, in the 3D subspace .
In the interval the invariant manifolds do not have two different forms. Along both directions the asymptotic curves have the same morphology. The asymptotic curves resemble the asymptotic curves of part B of the unstable invariant manifolds for (the case that we described in the section 4.3).

4.4 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 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 .
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 ) 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 . 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 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 x1v1 and x1v1 4D stability islands, which we can always see best in the 3D projection and the 2D projection . In all other projections the sticky regions are discernible as denser regions of the phase space embedded in clouds. The 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 transitions, as they are not expected to exist away from such a transition point or after a transition.
5 Double Instability
Double instability in rotating galactic potentials usually appears for large energies. In the x1 orbital tree (Skokos et al 2002) 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 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 zaxis orbits. These are p.o. along the rotation axis z in rotating triaxial potentials (Heisler et al. 1982).
5.1 x1v2
The periodic orbits of the x1v2 family are DU for large energies (Fig. 1). At first we study the case for . The periodic orbit has initial conditions . By adding even a tiny perturbation (e.g. or) to its initial conditions, along any direction, we find clouds of points in all 3D projections of the 4D space of section These are 4D clouds of points since we find scattered points in all 3D projections and mixing of colors when we apply the method of color and rotation. Their consequents occupy a volume in the 4D space of section, lying roughly in the region .
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 , for the same fixed value of the perturbation, . We find again 4D clouds of points as in the case of the double unstable periodic orbit at .
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 (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 and . The corresponding eigenvalues are and . The initial conditions , that we use for the computation of the asymptotic curves, are given by:
(4) 
We computed the asymptotic curves, in two cases which we named “I” and
“II” respectively. In the first case, we kept constant and equal to
zero and we used as a parameter. In the second case, we kept
and we varied . In case “I” we have taken values and
with a step
(2000 initial conditions). In the second case “II” we took values and
we varied in the interval
with a step (13500 initial conditions). In both cases “I” and
“II” we computed only 3 consequents for every initial
condition
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 4dimensional 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 and small 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 periodic orbits of the x1 tree. However, it is not necessarily typical for the dynamics in the neighborhood of 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 zaxis p.o.’s, which differs in many aspects from what we have found in the neighborhood of the x1v2 p.o.’s we presented up to now.
5.2 zaxis
As already mentioned, the orbits of the zaxis family lie entirely on the rotational, z axis, of the system. Thus, we consider now the surface as surface of section with . Earlier studies have shown that this family becomes double unstable in slowly rotating triaxal systems (Martinet & de Zeeuw 1988, Patsis & Zachilas 1990). We investigated the stability of this family in the potential (2) rotating very slowly with . For this value of the zaxis family has large regions of double instability (DU). In this case, as we vary the energy, the zaxis family for large negative values of is stable and for 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 yzplane. Then the zaxis family, for , becomes double unstable and it remains DU for larger values of .
In the 4D space of section the initial conditions of the zaxis orbits are . For we present the case of the orbit in the neighborhood of the zaxis p.o., which deviates from it in the direction by . We call this orbit “duz1”.
In the 2D 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 “”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 “”figure structure the symbol indicates the zaxis 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 projection is as we see in Fig. 19b.


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 “”figure structure of Fig. 19a are not on a smooth 4D surface. The application of the colorrotation method shows that we can hardly observe a smooth color variation along any direction. These are indications of a chaotic behavior.
The finite time index points also to the chaotic character of the orbit. For most of the time when the consequents in the (or the ) projection stay on the “”shaped figure, the index varies around 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 evolves as the consequents start filling all available phase space. It increases continuously up to and beyond.
The range of perturbations for which we find the “”figure is registered in Table 2. Outside the range given for the coordinate we find orbits sticky to rotational tori around the sao and its symmetric orbit. For absolutely larger perturbations in the , 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 zaxis 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 and . The corresponding eigenvalues are and . The initial conditions , which we use for the computation of the asymptotic curves, are given by:
(5) 
We calculate the manifold by using for the values with a step and for every value of we consider for the values again with a step . For every initial condition we compute 10 consequents. Then we have totally 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 “”figure in the 3D projection 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 “”shaped structure (Fig. 21a). In Fig. 21b we give the 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 “”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 and a step in (5), the asymptotic surface expands considerably also in the direction and its 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 “”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.

6 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 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 “”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 from the initial conditions of the simple unstable x1v2 for 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 at . We call this orbit “du1”.
In addition we include an orbit in the neighborhood of a zaxis p.o., i.e.

The orbit that builds an “”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 , we define as “diffusion time” . Then the “diffusion speed” is the ratio
The quantities and for the orbits in the neighborhood of a periodic orbit can be easily computed in two steps:

As time increases we compute the mean distance from the periodic orbit up to this time of all consequents and then we compute a mean volume . In other words we compute a “mean volume” of the hypersphere centered on the periodic orbit with radius the mean distance , i.e the quantity: .

After time the consequents occupy the maximum volume . This is represented by the maximum peak in a 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 for an interval , corresponding to 4200 intersections. We observe that for in our system units (indicated by the left arrow in Fig. 23a corresponding to 2580 intersections) 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 increases but the slope of the curve is smaller than the slope beyond the right arrow. The second arrow corresponds to time 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 projection (Fig. 5) sticking on the asymptotic curves of part B of the unstable manifold (section 4.4). Beyond the right arrow 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 corresponds to for . The speed of diffusion is . After the peak of in Fig. 23b the consequents arrive at distances from the periodic orbit that are smaller or equal to the maximum distance. This means that decreases a little and finally, it levels off to a constant value. This happens for all chaotic orbits.


The evolution of and the quantity 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 and . This means that the speed of diffusion is now , 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 for the orbit “du1” in the neighborhood of the p.o. x1v2 at . We observe that the curve has a maximum and . In this case the speed of diffusion is that is 3.65 times and 5133 times larger than the diffusion speeds of the orbits “u1” and “u” respectively. For “u1” and “du1” 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 zaxis p.o., i.e. for the orbit “duz1”, the evolution of the quantity is given in Fig. 25. Initially 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 “” structure, as in the 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 “”shaped structure. Thus, the mean volume decreases. Then the next consequents complete the “” structure by populating both lobes with similar number of consequents, so increases until , which corresponds to 3200 intersections (arrow numbered “2”). Then we observe a small plateau until the point indicated by the arrow labeled “3” (). The time at this point corresponds to 4200 consequents. This is the number of consequents that stay on the eightfigure surface before they depart from it and occupy a larger volume of the phase space. For larger times increases. In Fig. 25b we have a peak at for . The diffusion speed is .


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 variation.
7 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 (Skokos et al 2002) and zaxis (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 points around the p.o. with mixing of colors in their 4D representations. The two first kinds of orbital behavior appear close to the transition of the x1 family at the point of bifurcation of the x1v2 family and indicate the presence of the phenomenon of stickiness in the dynamics of the system.

The asymptotic curves of the x1v2 p.o. close to the transition of the central family of periodic orbits x1 are composed of two parts. Part A (calculated for in (3)) winds around deformed rotational tori belonging to orbits deviating from the x1 orbits in the or direction. Part B (calculated for 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 transition and a nearby 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 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 transition and in general there are no rotational tori associated with stable families in the neighborhood of the 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 zaxis 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 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 zaxis 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 zaxis 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.
Acknowledgments We thank Dr. C. Efstathiou for comments on an earlier draft of this paper.
8 References

Arnold, V.I. [1988], Geometrical Methods in the Theory of Ordinary Differential Equations (2nd edition, SpringerVerlag NY).

Broucke, R. [1969] “Periodic Orbits in the Elliptic Restricted ThreeBody Problem”, NASA Tech. Rep. 321360, 1125.

Contopoulos, G. [2002] Order and Chaos in Dynamical Astronomy (SpringerVerlag, New York Berlin Heidelberg).

Contopoulos, G. [2009] “Ordered and Chaotic Orbits in Spiral Galaxies” in ‘Chaos in Astronomy , G. Contopoulos and P.A. Patsis (eds) (SpringerVerlag, Berlin, Heidelberg), p.322.

Contopoulos, G. & Harsoula, M. [2008] “Stickiness in Chaos”, Int. J. Bif. Chaos 18, 29292949.

Contopoulos, G. & Harsoula, M. [2010] “Stickiness effects in chaos”, Celest. Mech. Dyn. Astron. 107, 7792.

Contopoulos, G. & Magnenat, P. [1985] “Simple threedimensional periodic orbits in a galactictype potential” Celest. Mech. 37, 387414.

Hadjidemetriou, J. [1975] “The stability of periodic orbits in the threebody problem”, Celest. Mech. 12, 255276.

Heisler, J., Merritt, D. & Schwarzschild, M. [1982] “Retrograde closed orbits in a rotating triaxial potential”, ApJ 258, 490498.

Katsanikas, M. & Patsis, P.A. [2011] “The structure of invariant tori in a 3D galactic potential”, Int. J. Bif. Chaos 21, 467496.

Katsanikas, M., Patsis, P.A. & Contopoulos, G. [2011] “The structure and evolution of confined tori near a Hamiltonian Hopf bifurcation”, Int. J. Bif. Chaos 21, 23212330.

Magnenat, P. [1981] Études numériques sur le computement des orbites stellaires dans des modeles galactiques a 2 ou 3 degrés de liberté PhD Thesis, Observatoire de Geneve.

Magnenat, P. [1982] “Numerical study of periodic orbit properties in a dynamical system with three degrees of freedom”, Celest. Mech. 28, 319343.

Martinet, L. & de Zeeuw, T. [1988] “Orbital stability in rotating triaxial stellar systems”, Astron. Astrophys. 206, 269278.

Miyamoto, M. & Nagai, R. [1975] “Threedimensional models for the distribution of mass in galaxies”, Publ. Astron. Soc. Japan 27, 533543.

Patsis, P. A. & Zachilas, L. [1990] “Complex instability of simple periodic orbits in a realistic twocomponent galactic potential”, Astron. Astrophys. 227, 3748.

Patsis, P.A. & Zachilas, L. [1994] “Using Color and rotation for visualizing fourdimensional Poincaré crosssections: with applications to the orbital behavior of a threedimensional Hamiltonian system”, Int. J. Bif. Chaos 4, 13991424.

Patsis, P.A., Skokos, Ch. & Athanassoula, E. [2002] “Orbital dynamics of threedimensional barsIII. Boxy/peanut edgeon profiles”, Mon Not. R. Astr. Soc. 337, 578596.

Skokos, Ch. [2010] “The Lyapunov Characteristic Exponents and their Computation”, Lect. Not. Phys. 790, 63135.

Skokos, Ch. [2001] “On the stability of periodic orbits of high dimensional autonomous Hamiltonian systems”, Physica D 159, 155179.

Skokos, Ch., Patsis, P.A. & Athanassoula, E. [2002] “Orbital dynamics of threedimensional barsI. The backbone of threedimensional bars. A fiducial case”, Mon. Not. R. Astr. Soc. 333, 847860.

Vrahatis, M.N., Isliker, H. & Bountis, T.C. [1997] “Structure and breakdown of invariant tori in a 4D mapping model of accelerator dynamics”, Int. J. Bif. Chaos. 7, 27072722.
Footnotes
 During the computation of asymptotic curves, we connect gaps that appear along them by considering more initial conditions in these gap regions.