Global dynamics in a self--consistent model of elliptical galaxy

In the present paper we study the global dynamics corresponding to a realistic model of self-consistent triaxial galactic system. We extend a previous work [17] where the authors investigate 3,472 orbits in this model at different energy levels, using Lyapunov exponents to measure chaoticity and frequency analysis to classify regular orbits. Here we first display the main properties of that potential and then focus our attention on the global dynamical features of the box domain for nine energy surfaces. Using the MEGNO as a fast dynamical indicator, we gain insight in the resonance structure at different energy levels, the way in which relatively large chaotic domains arise due to overlapping as well as crossings of resonances and we measure the fraction of chaotic motion in the energy space. It is interesting to notice that the flatness of the model varies over a rather wide range, namely from ~ $0.5$ to ~ $1$, and the fraction of chaotic motion ranges from ~ $0.15$ at small energies up to ~ $0.75$ at moderate values of the energy, decreasing then again down to values close to ~ $0.4$ where the system becomes nearly spherical.

Abstract. In the present paper we study the global dynamics corresponding to a realistic model of self-consistent triaxial galactic system. We extend a previous work [17] where the authors investigate 3,472 orbits in this model at different energy levels, using Lyapunov exponents to measure chaoticity and frequency analysis to classify regular orbits. Here we first display the main properties of that potential and then focus our attention on the global dynamical features of the box domain for nine energy surfaces. Using the MEGNO as a fast dynamical indicator, we gain insight in the resonance structure at different energy levels, the way in which relatively large chaotic domains arise due to overlapping as well as crossings of resonances and we measure the fraction of chaotic motion in the energy space. It is interesting to notice that the flatness of the model varies over a rather wide range, namely from ∼ 0.5 to ∼ 1, and the fraction of chaotic motion ranges from ∼ 0.15 at small energies up to ∼ 0.75 at moderate values of the energy, decreasing then again down to values close to ∼ 0.4 where the system becomes nearly spherical.
1. Introduction. One of the aims of galactic dynamics is to build self-consistent models of stellar systems, i.e., to get a mass distribution that produces a gravitational potential such that the orbits that it supports yield that same mass distribution. In simple cases, like spherical systems or disks, one can choose a certain potential, use the Jeans theorem (or its equivalent, the Boltzmann equation) in order to select a mass distribution compatible with that potential and, finally, with the Poisson equation make sure that such mass distribution creates the potential initially chosen. Analytical solutions are available, for example, for the Schuster (or Plummer) sphere and for the Mestel disk, although other models, like the King spheres, only admit numerical solutions (see, e.g. [2]). Self-consistency is much more difficult to attain in the triaxial systems needed to represent elliptical galaxies (see, e.g. [15]), and numerical methods are unavoidable to investigate them.
It is relatively simple to obtain models that resemble elliptical galaxies creating a random spherical distribution of point masses, with a very small initial kinetic energy, and following its collapse using an N-body code: one ends up with a system in virial equilibrium, with a de Vaucouleurs mass distribution, the more triaxial the smaller the initial kinetic energy (see [1]). Now, the number of stars in a real elliptical galaxy (say, 10 12 ) is much larger than the number of bodies that can be included in a numerical experiment (say, 10 6 ). Therefore, while the potential of the galaxy can be regarded as smooth and constant with time, the N-body potential is rather grainy, because the masses of the bodies are not negligibly small compared to the total mass of the system, and it changes slightly as time goes by, due to the more or less random changes in the N-body distribution as the bodies move along their orbits. Thus, the orbits that one would obtain in the N-body potential would be fairly different from the orbits in the galaxy that one is trying to model (e.g., the lack of smoothness and time constancy would yield spuriously large chaoticity).
The solution to this problem is to select a time after virial equilibrium has been reached and to fit the potential of the corresponding N-body distribution with an interpolating equation: the potential obtained from that equation will of course be smooth and constant, so that the orbits computed in it will mimic reasonably well those in a real galaxy. Moreover, since the selected N-body distribution is consistent with that potential, if one takes a randomly selected sample of those bodies and uses their positions and velocities as initial conditions to compute the orbits, then the different classes of computed orbits will appear in the same proportions as they are present in the galaxy. This method was pioneered by [23], who used it to investigate a rotating bar, and has been employed to model not only elliptical galaxies (see, e.g., [14] or [25]), but also galactic satellites [18].
In [17] the authors followed the collapse of a system of 10 5 bodies with a quadrupolar N-body code; after eliminating the bodies with positive and very low negative energies and letting the system relax again, they ended up with a system of 86,818 bodies resembling an elliptical galaxy (see their Fig. 1): their system follows a de Vaucouleurs law (see their Fig. 2) and it is strongly triaxial with flattening increasing from the border of the system to its center (see their Table I). They represented the potential with an equation fitted to an equilibrium N-body distribution and randomly selected a sample of 3,472 bodies to investigate the spatial distribution of the different classes of orbits in that potential. More recently, in [16] the author found that the system is not completely stationary, but that it experiences an extremely slow figure rotation (i.e., it rotates even though its total angular momentum is zero); the rotation is so slow (less than one revolution in a Hubble time) that no significant differences were found between the relative proportions of the different classes of regular orbits in the non-rotating and rotating systems. Nevertheless, probably because despite being extremely slow the rotation produces a symmetry break, there is a small but highly significant increase of 4.0% ± 1.2% of the chaotic orbits in the rotating system. Interestingly enough, Statler et al. [24] found some evidence of very slow figure rotation in NGC 4365 and found seemingly safe to conclude that it was dynamically unimportant in the observable part of that galaxy.
The potential introduced in [17] is thus very interesting, since it corresponds to a self-consistent N-body system that reproduces several features (such as mass distribution, flattening, triaxiality, figure rotation) of real galaxies. Nevertheless, no detailed studies of the properties of that potential has been carried out, so that in the present paper we undertake such an investigation. We deal here with orbits that go through the very center of the system; the investigation of other orbits (e.g., tubes), a detailed study of the dynamical consequences of the slow figure rotation, as well as the efficiency of the chaotic diffusion will be the subject of subsequent investigations.
In Section 2 we provide a brief summary concerning the MEGNO, while Section 3 is devoted to introducing the galactic model under study. In Section 4 we present the global dynamical properties of this system, for a particular but representative set of initial conditions.
2. The mean exponential growth factor of nearby orbits (MEGNO). Here we summarize the main features of the Mean Exponential Growth factor of Nearby Orbits (MEGNO hereafter) that is described in detail in [5]. This alternative tool to explore the phase space belongs to the class of the so-called fast indicators of dynamics.
Let H(p, q) with p, q ∈ R N denote an N -dimensional Hamiltonian, assumed to be autonomous just for the sake of simplicity since this is actually not required for the present formulation. On introducing the notation x = (p, q) ∈ R 2N , v = (−∂H/∂q, ∂H/∂p) ∈ R 2N , the equations of motion reaḋ Let γ(t) be an arc of an orbit of flow (1) on a compact energy surface and the full positive orbit is γ = lim t→∞ γ(t).
Relevant information about the flow in the vicinity of any orbit γ is gained through its largest Lyapunov Characteristic Number (LCN hereafter), σ(γ), defined as with δ(γ(t)) and δ 0 "infinitesimal displacements" from γ at times t and 0, respectively, (see below) and where · is some norm. The fact that the LCN measures the "mean exponential rate of divergence of nearby orbits", is stated explicitly when recasting (3) in the integral form with δ ≡ δ ,δ ≡ dδ/dt =δ · δ/ δ , the bar denoting time-average. Recall that the tangent vector δ satisfies the variational equatioṅ where Λ = Dv is the Jacobian matrix of the vector field v. We now introduce the MEGNO, Y (γ(t)), through the expression which is related to the integral appearing in (4). Notice that in the case of an exponential increase of δ, δ(γ(t)) = δ 0 exp(λt), the quantity Y (γ(t)) can be considered as a weighted variant of the integral in (4). Instead of the instantaneous rate of increase, λ, we average the logarithm of the growth factor, ln(δ(γ(t))/δ 0 ) = λt.
Let us now describe the MEGNO's asymptotic behavior in order to show how it serves to give a clear indication of the character of different types of orbits. In the first place, let us consider orbits on irrational tori for a non-isochronous system. As shown in [5], for any such (quasi-periodic) orbit, γ q , the temporal evolution of Y (γ q (t)) is given by where λ q is the linear rate of divergence around γ q and O denotes an oscillating term with zero average. The lim t→∞ Y (γ q (t)) does not exist but, on introducing the time-average it can readily be shown that Therefore, for the case of quasi-periodic motion, Y (γ) is a fixed constant, independent of γ.
For an irregular orbit, γ i , within any stochastic component, for which whereÕ is some oscillating term of bounded amplitude (which in general is neither periodic nor quasi-periodic, but has zero average). On averaging over an interval large enough we get Then, for a chaotic orbit both Y (γ i (t)) and Y (γ i (t)) grow linearly with time at a rate equal to the LCN of the orbit or one half of it, respectively. Only when the phase space has an hyperbolic structure does Y grow with time. Otherwise, it saturates to a constant value, even in the degenerated case in which δ grows with some power of t, say n, where Y → 2n as t → ∞. Therefore, the MEGNO's temporal evolution allows being summed up in a single expression valid for any kind of motion. The asymptotic behavior of Y (γ(t)) may be written in the fashion Y (γ(t)) ≈ a γ t + d γ , where a γ = σ γ /2 and d γ ≈ 0 for irregular, stochastic motion, while a γ = 0 and d γ ≈ 2 for stable quasi-periodic motion. Departures from the value d γ ≈ 2 indicate that γ is close to some periodic orbit, being d γ 2 and d γ 2 for stable and unstable periodic orbits, respectively (see [5] for details). Note also that the quantityσ 1 = Y /t verifieŝ (for quasi-periodic and chaotic orbits, respectively) which show that in the case of regular motionσ 1 converges to 0 faster than σ 1 (which goes to zero as ln t/t), while for chaotic motion both magnitudes approach the positive LCN at a rather similar rate. Further details on the MEGNO's performance when applied to the study of global dynamics in 2D Hamiltonians as well as the advantages of deriving the LCN from a least squares fit on Y , as well as a generalization of the MEGNO and its application to discrete dynamical systems are given in [4], [5] and [7]. An interesting application to a 2.5D problem, namely, the classical Arnol'd's diffusion problem when the two small parameters are equal, may be found in [21]. The associated splitting of separatrices has been studied in [22]. Recent applications of the MEGNO to Solar System dynamics may be found in, for instance [3], [5], [9], [10], [11] and [12].
3. The model. In the present work we study the potential given in [17] Φ being p the softened radius given by and by In any case, the softening parameter, ǫ, is equal to 0.01. The functions f s (p) were computed from an N-body simulation and fitted with equations of the form (see [17]): The adopted values for C s , k s , q s , l s are the ones presented in Table 1.  Table 1. Values for the coefficients of the functions f s given by Eq. (16).
The stationary character for the values of the parameters given in Table 1 has been tested by performing several fits at different times after virialization, resulting with a precision of the order of 0.1%. Fig. 1 displays the behaviour of the f s 's with r, being f z < 0, f 0 , f x > 0, and f x > |f z | for all values of r while for r 0.36 it is f 0 > f x . There we plot the functions f s with the concomitant sign with which they appear in Eq. (13). The resulting potential is triaxial with semi-axes X, Y, Z satisfying the condition X > Y > Z. The minimum of the potential, which is very close to −7, takes place at the origin. The potential is, of course, less flattened than the mass distribution (see Table I of [17]).
The semi-axis ratios of the potential are represented in the upper-left panel of Fig. 2. In any case the ratios lie in the range 0.5 < Y /X, Z/X ≤ 1. For energy values within the interval −1 < E < 0, both ratios approach each other, and for E close to 0 they become rather similar and close to 1. This is the expected behaviour, since in this energy range r could reach large values and, consequently, both f x and f z become rather small, so that the system approaches spherical symmetry (see below). The minimum value for a semi-axis ratio appears close to E = −6 where Z/X ≈ 0.52. The maximum distance to the center at each energy surface, i.e. the long axis X, is represented in the upper-right panel of Fig. 2, where it can be seen that for large energies, around E −0.5, r could reach values larger than 2. Although the triaxiality is usually defined as:  Figure 1.
Dependence of the f s given by Eq. (16) with the parameters given in Table 1 on r.  Figure 2. Semiaxis ratios Y /X, Z/X and long axis X as functions of the energy (left and right top panels respectively). Triaxial parameter T , and T * given by (17) vs. the energy (left and right bottom panels). Fig. 2) we prefer to adopt as a measure of the triaxiality of the potential the magnitude which is represented at the right-bottom of Fig. 2. This magnitude is well defined when Z → Y → X, that is when E → 0, which is not the case for T . Let us notice that along the energy range, T * has an increasing behaviour, reaching its maximum at E ≈ −3, and then decreases to 0 for E → 0, as expected.
In terms of spherical coordinates the potential, can be recast as: where, withf s (r) = f s (p(r)); s = 0, x, z. The dependence of the coefficients with r is shown in Fig. 3. Compared to V 0 , the other V s 's are seen to be significant only for r 0.3. Thus, considering the smallness of the coefficients V s , s = 1, 4, the Hamiltonian can be written as a near-integrable one in the form: V (r) = V 1 (r) cos 2φ + V 2 (r) cos 2θ + V 3 (r) cos 2(θ + φ) + V 4 (r) cos 2(θ − φ), (22) where p r =ṙ, p θ = r 2θ , p φ = r 2φ sin 2 θ.
H 0 is an integrable Hamiltonian being the three unperturbed action-like integrals, whileV could be considered as a small perturbation. Therefore, when we refer to unperturbed motion, we mean orbits in H 0 . The unperturbed system is just a central field and, thus, the motion is confined to a plane. Let us take θ = π/2, which means the xy plane. Therefore, p θ = 0 and I o 2 1 = I o 2 , and r oscillates between two boundaries, r m (E 0 , I o 1 ) ≤ r o (t) ≤ r M (E 0 , I o 1 ), with frequency ω r , while φ o varies onto the circle S 1 . The frequency in the tangential direction is ω φ = κ ω r where κ = ∆φ/2π < 1 is, in general, irrational (see [2]). The time evolution of φ o can be written as φ o (t) = ω φ t+ϕ(t) where ϕ is a 2π/ω r -periodic function of time.
Let us simplify the problem restricting the motion to the invariant planes xy, xz and yz. For the xy plane we set, as before, θ = π/2, p θ = 0 and the system of Eqs. (27), (28) reduces to the single ordinary differential equation: dp φ dφ p φ = 4r 2 V 1 (r) sin 2φ.
Following [7], the latter can be written as: where g and g 1 (r(t)) are the average and oscillating parts of 4r 2 (t)V 1 (r(t)) respectively: g = 4V 1 (r)r 2 < 0, g 1 (r(t)) = 4V 1 (r(t)) r 2 (t) − g, (see [7]). To keep up to the first order of the perturbation, in the second term in the first of (30) we can replace the actual values of r(t), φ(t) by the concomitant unperturbed values r o (t), φ o (t). Since the unperturbed motion is completely regular, we can expand g 1 (r o (t)) sin 2φ o (t) in Fourier series, with fundamental frequencies ω r and ω φ whereg k are certain complex coefficients. Assuming quasi-periodicity (which is the prevalent behaviour if the perturbation is small), i.e. κ irrational, we easily see from (32) that g 1 sin 2φ o ≈ 0. Hence, if we average the first equation in (30) over several radial periods we see that is an approximate invariant.K plays the role of the total energy in a simple pendulum model where Ω 1 is the small oscillation frequency. Therefore, there are two critical values ofK, namely, −Ω 2 1 and Ω 2 1 . ForK = −Ω 2 1 , (φ, p φ ) = (0, 0) we have a stable equilibrium point: the motion is stable along the x axis. On the other hand, forK = Ω 2 1 we have the separatrix and the unstable equilibrium points are (φ, p φ ) = (±π/2, 0); the motion along the y axis is thus unstable. In the xy plane, the domain of box orbits that oscillate about the long-axis corresponds to |K| < Ω 2 1 and the domain of loop orbits, that rotate about the origin, toK > Ω 2 1 . The separatrix, p s φ = ±2Ω 1 cos φ s , separates then different kinds of motion: oscillations and rotations; i.e. box and loop orbits. ForK ≫ Ω 2 1 , it isK ≈ p 2 φ /2: the kinetic energy in the tangential direction. The largest value ofK corresponds to the largest p φ , which appears for the 1:1 periodic orbit in the xy plane. For V 1 = 0 but small, this periodic orbit should not be circular but elliptical with a small eccentricity.
Since Ω 1 is a measure of the amplitude of the torque, we conclude that the relevant parameter that defines the orbital family in the plane is the relative value of the rotational energy with respect to the strength of the torque, which in turn depends on how elongated the potential is.
From the above discussion it turns out that a limit angle, φ l , could exist, which is another way to conclude that |K| < Ω 2 1 for boxes. However it is important to remark that this bound for φ appears for r bounded away from 0. When p φ is small, which is the case for boxes, the analysis of the motion in a neighbourhood R of r = 0 should be done in a different fashion since the origin is a singular point in this description. Since the potential is regular at r = 0, we can approximate Φ(x, y, 0) by a harmonic oscillator in R. The approximate invariants of motion are then the energy in each degree of freedom E x , E y . But the Lissajous-like orbits in a harmonic oscillator with incommensurable frequencies are dense in R whenever E x , E y = 0, so no bound for φ exists while the star is in R.
Thus, taking into account Eq. (20) we obtain the following approximate invariants: whereḡ is the average part of for Eq. (34), and of r 2 (t)(V 2 (r(t)) − 2V 3 (r(t)) for Eq. (35). Considering the sign in the second term in Eq. (36), (θ, p θ ) = (π/2, 0) is a stable equilibrium point: the motion is stable along the x axis. On the other hand, forĴ = Ω 2 2 the unstable equilibrium points are (θ, p θ ) = (0, 0): the motion along the z axis is unstable. The domain of box orbits, that oscillate about the longaxis in the xz plane, corresponds to |Ĵ| < Ω 2 2 and the domain of loop orbits, that rotate about the origin, toĴ > Ω 2 2 . For φ = π 2 and p φ = 0 the motion is restricted  to the invariant plane yz and the analysis is the same as that given above. Thus (θ, p θ ) = (π/2, 0) is the stable equilibrium point: the motion is stable along the y axis and the unstable equilibrium points are (θ, p θ ) = (0, 0): the motion along the z axis is unstable.
In Fig. 4 we present the stability analysis of the three axial periodic orbits by means of the MEGNO, Y , (see [5], [7], [8]). The value of the MEGNO for stable orbits should approach asymptotically to the value 2, while in the case of unstable orbits, its value grows linearly with time, thus its final value depends on the final time and should be very large. However, in this case we observe a small window of stability for the x and y periodic orbits for E −5.5 where the MEGNO is close to 0. This is the expected value of Y for stable motion in nearly linear systems. Since for small energy values the potential should behave as a perturbed harmonic oscillator, the computed values of the MEGNO for rectilinear periodic motion indicate stability for these two orbits. Notice, however, that the x-axis periodic orbit becomes again stable for E −0.3. On the other hand, the z-axis periodic orbit is always unstable, fully consistent with the above analysis of stability. From these figures we could infer the energy values for which the tube family arise, from bifurcation of the x and y periodic orbits. These bifurcation occurs very close to E = −6 and to E −5.
Out of intuition, we may think that box orbits respect the global invariant H 0 , and the local onesK andĴ. The full problem in 3D seems to be rather difficult to be treated analytically. In the same direction, one may understand the existence of three families of tube orbits, just by analysing the stability of the 1:1 unperturbed circular orbit in each invariant plane (see for instance, [19]). 4. Global dynamics on several energy manifolds. In the previous work [17] the authors studied 3, 472 orbits in different energy surfaces whose distribution is shown in Fig. 5. There, the authors perform an orbital classification (in boxes and tubes) and also investigate chaoticity by means of the Lyapunov Characteristic Numbers (LCN) in order to discriminate between what they call fully chaotic, partially chaotic or regular orbits. By fully chaotic orbits they mean those having two positive LCNs, and by partially chaotic ones, those that present only one positive LCN.
In the present work we investigate the global dynamics onto several energy surfaces using the MEGNO as a fast indicator. In each case we have considered about 52, 000 orbits, except for large energies, for which more than 120, 000 orbits have been analysed. We have adopted a total motion time of 3, 000 units and integrated the equations of motion together with their first variationals, the initial conditions for the variationals being taken at random in phase space and with unit norm. We have used a Runge-Kutta 7/8 th order integrator (the so-called DOPRI8 routine, see [13]). The accuracy in the preservation of the value of the energy is ∼ 10 −13 . In this first paper we address the analysis of the box domain, this family of orbits having a natural space onto which be represented. Taking as initial conditions in configuration space x = y = z = 0, and since the Hamiltonian is a global invariant, an adequate space to investigate the global dynamics on this submanifold of phase space is, for instance, the plane p 2 x p 2 z . In a forthcoming paper we will be dealing with the tube domain, families for which the global dynamical properties should be investigated in different spaces (see the classical paper [20]).
We deal with the normalized values of p 2 s obtained dividing them by the maximum value they can reach at each energy level. We present the obtained values for the MEGNO in a contour-like plot for different energy levels, namely, for E = −6, −5, −4, −2 in Fig. 6 and for E = −1, −05, −0.3 and −0.1 in Fig. 7. (The plot corresponding to E = −3 has not been included since it is not substantially different from that for E = −2, which is in concordance with what can be observed in Fig. 8.) There, the values of the < Y > were binned in three intervals. Points depicted in white correspond to values of < Y > in the interval [0, 2.02), that is, to regular/stable dynamics; while those points depicted in green correspond to (local) mild unstable dynamics. Red points identify (local) strong unstable motion.
Each one of these plots let us gain insight in the resonance structure of the system for a particular value of the energy, evincing the arise of relatively large chaotic domains due to overlapping as well as crossings of resonances.  Figure 8. Fraction of chaotic orbits on the energy space obtained using the MEGNO to separate regular from chaotic orbits.
Let us turn our attention to Fig. 6. Notice that for a low value of the energy, such as E = −6, the regular regime clearly prevails and there we observe a strip of chaotic motion near p 2 x + p 2 z = 1. A wider chaotic domain shows up for E = −5, where the rather intricate crossing of resonances can be clearly visualized, the main chaotic regime lying close to the p 2 z axis. The chaotic domain is seen to increase as larger values of the energy are considered, reaching its maximum widespread for E = −2. Fig. 7 displays the concomitant plots for large values of the energy. Intersections of resonances can clearly be distinguished. Note that for E = −0.5 several small stability regions are embedded in the mild chaotic component. As the energy is increased, the chaotic domain shrinks, becoming rather compact and weaker for E = −0.1, where regular motion prevails.
In Fig. 8 we present the fraction of chaotic motion as a function of the energy, adopting a threshold of 2 for < Y >, below which orbits are regarded as regular (similar results were obtained using the threshold 2.02). Note that there is some value of E, close to −3 where a transition to nearly global stochasticity takes place, in the sense that chaotic motion prevails in phase space, the maximum occurring for E = −2 where nearly 75% of the energy surface is covered by chaotic motion. It is clear from the figure that this effect reverses for E > −0.5, where the regular component is the one that predominates as larger energy levels are considered. It should be mentioned that the fraction of chaotic motion over each energy surface has been computed considering about 52, 000 orbits for E ≤ −1 and nearly 120, 000 orbits for E > −1, except for E = −6.5 where only a few thousand orbits have been taken into account.

5.
Discussion. Let us turn back to Fig. 3, where the magnitude of the different coefficients appearing in the multipolar expansion of the potential is displayed. There we observe that the non-integrable perturbation to the central field model is significant only for 0 < r 0.3. On the other hand, the right upper panel in Fig. 2 shows that the values r ≤ X 0.3 correspond to energies in the range −4 E −2. Therefore, we might expect that the strongest effect of the perturbation occurs for such energy range, since for those values of the energy any particle will be interacting most of the time with a highly non-integrable field. This is in fully agreement with the results presented in the previous section, namely in Figs. 6, 7 and 8.
In terms of macroscopical parameters, the maximum value of the triaxiality measured by T * (right bottom panel in Fig. 2) is reached around E ∼ −3. Thus we note that the parameter T * provides a good indication of the strength of the perturbation at different energy levels. Indeed, as Fig. 6 shows, for −4 E −2 the energy surface is dominated by chaotic motion mainly due to the overlap of resonances.
For small energies, E −6, the system could be approximated by a perturbed harmonic oscillator, the small unstable, chaotic regions that appear at these energy levels arising from the unstable character of the z-axis periodic orbit. In the range −6 E −5 both the x and y-axis periodic orbits become unstable and larger domains of chaotic motion should appear, as Fig. 6 reveals. The small fraction of the energy surface covered by chaotic motion at large energies is certainly due to the fact that the system becomes more spherical and the x-axis periodic orbit approaches stability leading to a larger region of regular motion.
From Fig. 8 we observe that the fraction of chaotic motion on the energy surface has a smooth behaviour, reaching values close to 0.75 for the E ≈ −2. This fraction decreases to 0.4 for E ≈ 0, and to 0.15 for E = −6.5.
The global dynamical properties of other starting spaces (e.g., tubes orbits) will be addressed in a forthcoming paper. Anyway, having attained the global dynamical structure corresponding to box orbit domain (which are the dominant orbits in elliptical galaxies), allows us to proceed, in a subsequent investigation, with the study of chaotic diffusion on phase space. In this direction, new interesting results concerning diffusion on phase space in conservative dynamical systems are presented in [6].