A&A 455, 499-507 (2006)
DOI: 10.1051/0004-6361:20054721

Global dynamics in galactic triaxial systems. I

P. M. Cincotta^1,2 - C. M. Giordano^1,2 - M. J. Pérez^1,2

1 - Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque S/N, La Plata (1900), Argentina
2 - Instituto de Astrofísica de La Plata (IALP), Consejo Nacional de Investigaciones Cientícas y Técnicas de la República Argentina

Received 19 December 2005 / Accepted 18 April 2006

Abstract
We present a theoretical analysis of the global dynamics in a triaxial galactic system using a 3D integrable Hamiltonian as a simple representation. We include a thorough discussion on the effect of adding a generic non-integrable perturbation to the global dynamics of the system. We adopt the triaxial Stäckel Hamiltonian as the integrable model and compute its resonance structure in order to understand its global dynamics when a perturbation is introduced. We provide a theoretical discussion about diffusive processes taking place in phase space.

Key words: methods: numerical - galaxies: structure - chaos - diffusion - galaxies: kinematics and dynamics

   
1 Introduction

Observations with HST revealed the presence of very high stellar densities at the centers of early-type galaxies (Crane et al. 1993; Ferrarese et al. 1994), suggesting a power law ( $r^{-\gamma}$) fit for them. The evidence for large central masses was also reforced from high-resolution kinematic studies of nuclear stars and gas, which revealed the presence of compact dark objects with masses in the range of $10^{6.5}{-}10^{9.5}~M_{\odot}$, presumably super-massive black holes (Ford et al. 1998). These observational results have produced a substantial change in the classic thinking on dynamics in triaxial galaxies.

Results obtained from numerical simulations show that the addition of a central mass to an integrable triaxial potential has strong effects on its dynamics, at least for the boxlike orbits that cover the central region of triaxial galaxies. Black holes and central density cusps scatter these particular orbits during each close passage inducing chaos in the system. The sensitivity of boxlike orbits to deflections also leads to a rounder central distribution of mass (see for instance Gerhard & Binney 1985, and Udry & Pfenniger 1988). This slow evolution towards axisymmetry suggests that stationary triaxial configurations cannot exist for a central density cusp. Valluri & Merritt (1998) find that in most of early-type galaxies, chaotic evolution is determined by the mass of the central black hole, ( $M_{\rm bh}$), rather than by the slope of the density profile. They show that when the central mass contains $2 \%$ of the galaxy mass, a transition to global stochasticity occurs. For such large values of  $M_{\rm bh}$, the box-orbit phase space is almost completely stochastic and diffusive processes could take place over very short timescales.

This result is substantially attractive because this critical black hole mass was close to the observed one (Kormendy & Ritchstone 1995) and also close to the mass that induced a sudden evolution toward axisymmetry in N-body simulations (Merritt & Quinlan 1998). However, this is no longer true, since from the works of, for instance, Ferrarese & Merritt (2000) and Merritt & Ferrarese (2001) it is known that the mass of black holes in galaxies, derived from the black hole demographic relationships, are $0.1{-}0.2\%$ of the ellipsoid mass in which they are embedded.

Merritt & Fridman (1996) arrive at the same conclusion analyzing two triaxial power law models: the steep ($\gamma =2$) and the weak ($\gamma =1$) cusp. They find, in agreement with Gerhard & Binney (1985) and Schwarzschild (1993), that triaxial galaxies with such high concentration of mass would evolve toward a central axisymmetry, as box orbits gradually lose their distinguishability. For these models, in which a large fraction of phase space is dominated by chaotic dynamics, the construction of self-consistent solutions requires the inclusion of stochastic orbits as well as the regular ones. A system thus built evolves, mainly close to its center, as stochastic orbits mix, through phase space. In order to obtain stationary solutions, they build the "fully-mixed models'' keeping stochastic orbits out of the central part of the system where chaotic orbits mix ergodically driving to a rounder distribution and destroying the triaxial self-consistency. They find that although it is possible to build this kind of solution for a weak cusp model, this is not the case for a strongly concentrated model. This would imply that triaxility is not consistent with a high central density. The discussion of the existance of stationary non-axisymetric stellar systems is still open. Poon & Merritt (2002) were able to construct triaxial equilibrium with central black holes which were both regular and stable, but under a very strong assumption. They need the chaotic building blocks to be fully mixed for the triaxial equilibria to co-exist with central singularities, and there is no evidence yet that in 3D systems with divided phase space, a completely connected chaotic component actually exists. Moreover, it seems that this could happen only when the chaotic component has a large measure ($\sim$1) and $t\to\infty$, which, from a physical point of view, would not be possible in galactic systems (see for instance Giordano & Cincotta 2004 and references therein and also Sect. 5).

All these conclusions rest on the strong hypothesis that diffusive processes drive stochastic orbits to mix, covering all chaotic regions of phase space on a timescale of the order of the Hubble time. Taking into account the fundamental role that diffusive processes should play in galactic dynamics (as well as in asteroidal and planetary dynamics), we will discuss the relevance of diffusive processes in phase space theoretically as well as numerically. In this first paper we address the theoretical discussion; numerical studies will be presented in a forthcoming paper.

Global dynamics of triaxial galactic systems have been studied in previous works, for example by Wachlin & Ferraz-Mello (1998) and Papaphilippou & Laskar (1998). They represented elliptical galaxies using, respectively, a generalization of a double-power-law spherical mass model and the axisymmetric softened version of the 3D logarithmic potential. For these models, they applied frequency map analysis (Laskar 1993) in order to investigate the orbital structure of the system. Both analyses revealed stochastic motion in both models. Although they show the main resonances of these systems and investigate the relevance of chaotic motion in these models, they did not focus on the study of the diffusive mechanisms that may present in these systems. Moreover, the frequency map analysis is not an efficient tool for such of studies since it does not provide a good measure of chaos (see Cincotta & Simó 2000; Cincotta et al. 2003).

Here, we compute the theoretical resonant structure of an integrable triaxial system. Particularly, we focus our attention on the Stäckel triaxial model, which we assume to be a rough approximation of an elliptical galaxy. Then we add a generic non-integrable perturbation in order to show how the resonance structure is distorted by the effect of this perturbation.

In Sect. 2, we briefly recall the main characteristics of the integrable triaxial Stäckel model. Numerical details related to the computation of the resonances and the analysis of the resonance structure are presented in Sect. 3. In Sect. 4, we review and provide a theoretical discussion about the transition to chaos as a consequence of resonance interaction. Finally, we discuss the different processes that could lead to chaotic diffusion in phase space.

   
2 The integrable model

Motion in a smooth gravitational field becomes quite simple if the number of isolating integrals equals the number of degrees of freedom, and much work in galactic dynamics has focused on finding integrable models for galactic potentials. Kuzmin (1956, 1973) showed that there is a unique, ellipsoidally stratified mass model for which the corresponding potential has three global integrals of the motion, quadratic in the velocities. Kuzmin's model, explored in detail by de Zeeuw (1985) who called it the "Perfect Ellipsoid'', has a large, constant-density core in which the orbital motion is similar to that of a three-dimensional linear oscillator. Every orbit in the core of the Perfect Ellipsoid fills a region close to a rectangular parallelepiped, or box. These trajectories were called box orbits. At higher energies in the Perfect Ellipsoid, box orbits persist and three new orbital families appear: the tubes (inner and outer long axis tubes and short axis tubes), orbits that preserve the direction of their circulation around either the long or short axis of the figure. Tube orbits respect an integral of motion analogous to the angular momentum, and hence - unlike box orbits - avoid the center. As pointed out by Schwarzschild (1981), these four families, box, inner and outer long axis tubes and short axis tubes, reproduce the general form of triaxial galaxies. The Perfect Ellipsoid model has been used by several authors as a first approximation to represent an elliptical galaxy.

In this section, we address some relevant aspects of the integrable Stäckel model - separable in ellipsoidal coordinates - whose three global integrals of motion admit analytical expressions.

The Perfect Ellipsoid is represented by stratified density function distributed over concentric ellipsoids of semi-axes ma, mb and mc, given by

$$% \rho={{\rho_{0}}\over{(1+m^{2})^{2}}},$$ (1)

where $\rho_{0}$ represents the central density and m,

$$% m^{2}={{x^{2}}\over{a^{2}}}+{{y^{2}}\over{b^{2}}} +{{z^{2}}\over{c^{2}}}, \quad \quad a \geq b \geq c \geq 0,$$ (2)

is constant on an ellipsoidal shell. Following Chandrasekhar (1969), it is possible to obtain the potential for the Perfect Ellipsoid, which in ellipsoidal coordinates  $(\lambda, \mu, \nu)$ adopts the Stäckel form:

$$% V(\lambda, \mu, \nu)=-{{1}\over{4 h_{\lambda}^{2}}}{{G(\lam... ...eta)}} -{{1}\over{4 h_{\nu}^{2}}}{{G(\nu)}\over{(\nu+\beta)}},$$ (3)

where $\alpha=-a^{2},~ \beta=-b^{2},~ \gamma=-c^{2}$ and  $h_{\lambda}^{2},~ h_{\mu}^{2},~ h_{\nu}^{2}$ are the metric coefficients of the ellipsoidal coordinates and $G(\tau)$ is a function given in terms of an elliptical integral of the third kind (see de Zeeuw 1985).

The Stäckel model has three explicit global analytic integrals, namely the adelphic integrals I₂ and I₃ and the total energy H, and indeed all orbits in it can be determined in terms of simple quadratures.

The adelphic integrals can be considered as generalizations of the angular momentum integrals that exist in axisymmetric and spherical potentials, but also as generalizations of the energy integrals always present in separable potentials in Cartesian coordinates. Every orbit in Stäckel potential is the sum of three motions, one in each coordinate. As a result, motion is bounded by coordinate surfaces.

The adelphic integrals I₂ and I₃ are, in fact, linear combinations of other integrals J and K:

$$% I_{2}={{\alpha^{2} H+ \alpha J+ K}\over{\alpha-\gamma}}, \q... ...\quad I_{3}={{\gamma^{2} H+ \gamma J+ K}\over{\gamma-\alpha}},$$ (4)

where the energy H, J and K are functions of the ellipsoidal coordinates and conjugate momenta (see de Zeeuw 1985).

As a consequence of the separability of the Stäckel potential, the equations of motion can be recast as:

$$% p_{\tau}^{2}={{H-V_{{\rm eff}}(\tau)}\over{2(\tau+\beta)}},$$ (5)

where $\tau= \lambda, \mu, \nu$, and $V_{{\rm eff}}(\tau)$denotes the effective potential,

$$% V_{{\rm eff}}(\tau)={{I_{2}}\over{(\tau+\alpha)}} +{{I_{3}}\over{(\tau+\gamma)}}-G(\tau).$$ (6)

This characteristic of the system allows us to analyze the orbital structure according to the values of  $p^{2}_{\tau}$ for each set of chosen values of the integrals. Since the motion is only possible for  $p^{2}_{\tau}\ge 0$, $\tau= \lambda, \mu, \nu$, the motion in each coordinate is either an oscillation between turning points defined by $p^{2}_{\tau}=0$, or a rotation whenever  $p^{2}_{\tau}>0$ for all $\tau$. The combination of these two kinds of motion will determine the class of orbit. The ranges where  $p^{2}_{\tau}$ are non-negative define a volume space, determined by the integrals, where the motion takes place. Thus, the orbital structure of the system may well be studied by analyzing the values of (5) for different combinations of chosen values for its integrals. Such a systematic study can be performed by fixing one of the integrals to compute the momenta in the space determined by the two remaining ones.

de Zeeuw (1985) carried out a detailed classification of orbits for this Stäckel model by analyzing the tridimensional space of integrals. A projection of this space on the energy surface is shown in Fig. 1, where the two inner curves in the figure separate the four different orbital families: box, inner and outer long axis tubes and short axis tubes. These curves are separatices of different kinds of motion and they play a central role in the dynamics of the system when it is perturbed (see Sect. 4).

Although the properties of this model have already been thoroughly studied (de Zeeuw & Lynden-Bell 1985), we will complement them by computing and analyzing the resonance structure of the system in its integrals space. The knowledge of this resonance structure will allow us to predict the global dynamical behavior of the system when a non-integrable perturbation is introduced.

	Figure 1: Orbital Classification for the Perfect Ellipsoid in the I2I3-plane, for a fixed value of the energy E = -0.3, de Zweeu (1985). Curves drawn in this plane are separatrices of different families of orbits.
Open with DEXTER

   
3 Resonance structure of the integrable model

In order to understand how an integrable system is modified by the effect of a non-integrable perturbation, we require a picture of its resonance structure, since it yields a preliminary knowledge of the way in which resonances in the integrable system may be distorted when a perturbation is turned on, avoiding the integration of the equations of motion and the variational ones. Further, this gives us information about the mechanisms that could give rise to the transition from regularity to chaos by different processes.

In the present section, we compute the resonance structure for the integrable Stäckel model in the non-perturbed integral space I =(H, I₂, I₃). At this stage we introduce action-angle variables, which in many respects are the natural variables for the description of the motion. Indeed, from the geometrical point of view the actions (which are functions of the global integrals) define the tori structure that foliate the phase space. For given values of the action, the angle variables describe the orbit on the torus.

Thus, we briefly outline the procedure for the case of motion separable in ellipsoidal coordinates, whose results apply not only to the particular case of the Perfect Ellipsoid, but to every Stäckel potential. The action variables $J_{\tau}$ associated with each ellipsoidal coordinate are defined by

$$% J_{\tau}={{1}\over{2\pi}}\oint p_{\tau} {\rm d}\tau= {{2}\o... ...\rm max}}} p_{\tau} {\rm d}\tau, \qquad \tau=\lambda,\mu,\nu,$$ (7)

where, from Eqs. (5) and (6), the momenta $p_{\tau}$ can be written (de Zweeu 1985)

$$% p_{\tau}^{2}= {{H-I_{2}/(\tau+\alpha)-I_{3}/(\tau+\gamma)+G(\tau)}\over{2(\tau+\beta)}},$$ (8)

and the integration is performed over all values of $\tau$ for which ${p_{\tau}}^2 \geq 0$.

From Eqs. (7) and (8) we obtain $J_{\tau}=J_{\tau}(H,I_2,I_3)$and denoting the action vector by $\bm {J}=(J_{\lambda}, J_{\nu}, J_{\mu})$, we can formally invert this relation in order to get H(J). For given values of the integrals, H=E, I₂=i₂ and I₃=i₃, we get the values of the actions $J_{\tau}$ that fix the torus where the motion takes place. The associated angle variables  $\theta_{\tau}$, canonically conjugate to the actions $J_{\tau}$, are linear functions of the time t,

$$% \theta_{\tau}={\omega}_{\tau}(\bm{J})\cdot t +\theta_{\tau}(0),$$ (9)

where $\theta_{\tau}(0)$ are constants determined by the three remnant initial conditions and the quantities (integrals) $\omega_{\tau}(\bm{J})$ given by

$$% \bm{\omega}(\bm{J})={\partial H(\bm{J})\over\partial\bm{J}}$$ (10)

are the frequencies of the motion in each degree of freedom.

If the three frequencies are incommensurable, the orbit densely fills the volume allowed to it by the values of the integrals of motion, i.e. the region where  $p_{\tau}^2$ are non-negative. Thus, for these particular values of the actions (or the integrals) we have a strong non-resonant 3D torus on which the motion is ergodic. This occurs when the frequencies satisfy the so-called diophantine condition (see, for instance, Giorgili 1990) which states that for a given integer vector m the frequency vector satisfies

$$% \vert\bm{m}\cdot\bm{\omega}\vert\ge\gamma~\vert\bm{m}\vert^... ...ert m\vert=\vert m_{1}\vert+\vert m_{2}\vert+\vert m_{3}\vert,$$ (11)

where $\gamma>0$ and $\alpha>2$ in case of a 3D system. Whenever a relation of the form

$$% \bm{m}\cdot\bm{\omega}(\bm{J})=0, \qquad \bm{m}\in\mathbb{Z}^{3}/\bm{0}$$ (12)

is satisfied, we obtain a resonance condition for the actions or the frequencies. Those values of the actions that fulfill (12) lead to a resonant torus, that is, orbits are not ergodic on a 3D torus but on a 2-dimensional one. It may be interpreted that the resonance condition provides a relation between the actions that lead to a new local integral which confines the motion to a manifold of lower dimensionality.

In order to obtain the resonance structure on an energy surface, a particular value of the energy H has to be fixed and the resonance condition (12) solved for each value of J.

The frequencies for ellipsoidal coordinates are given by (10) and, since actions do not admit explicit analytical expressions in terms of the integrals, frequencies (10) have to be calculated by numerical means. Indeed, since $H=H(J_\lambda,J_\mu,J_\nu)$, we may write (de Zweeu 1985)

 

$$% {{\partial H}\over{\partial H}} \omega_{\lambda}{{\partial J_{\lambda}}\over{\partial H}}+ \omega... ..._{\mu}}\over{\partial H}}+ \omega_{\nu}{{\partial J_{\nu}}\over{\partial H}}=1,$$ =

$${{\partial H}\over{\partial I_{2}}} \omega_{\lambda}{{\partial J_{\lambda}}\over{\partial I_{2}}}+ \o... ...over{\partial I_{2}}}+ \omega_{\nu}{{\partial J_{\nu}}\over{\partial I_{2}}}=0,$$ = (13)

$${{\partial H}\over{\partial I_{3}}} \omega_{\lambda}{{\partial J_{\lambda}}\over{\partial I_{3}} }+ \... ...over{\partial I_{3}}}+ \omega_{\nu}{{\partial J_{\nu}}\over{\partial I_{3}}}=0.$$ =

Then,

 

$$% \omega_{\lambda}={{1}\over{\Delta}}~ {{\partial(J_{\mu},J_{\nu}... ...\over{\Delta}}~ {{\partial(J_{\lambda},J_{\mu})}\over{\partial (I_{2},I_{3})}},$$ (14)

where we have written

$$% \Delta= {{\partial(J_{\lambda},J_{\mu},J_{\nu})}\over{\partial (H,I_{2},I_{3})}}\cdot$$ (15)

The partial derivatives that occur in the above expressions can be evaluated by differentiating under the integrals, (7)

 

$$% {{\partial J_{\tau}}\over{\partial H}} {{2}\over{\pi}} \int_{\tau_{{\rm min}}}^{\tau_{{\rm max}}} {{{\rm d}\tau}\over{(\tau+\beta)p_{\tau}}},$$ =

$${{\partial J_{\tau}}\over{\partial I_{2}}} -{{2}\over{\pi}} \int_{\tau_{{\rm min}}}^{\tau_{{\rm max}}} {{{\rm d}\tau}\over{(\tau+\alpha)(\tau+\beta)p_{\tau}}},$$ = (16)

$${{\partial J_{\tau}}\over{\partial I_{3}}} -{{2}\over{\pi}} \int_{\tau_{{\rm min}}}^{\tau_{{\rm max}}} {{{\rm d}\tau}\over{(\tau+\beta)(\tau+\gamma)p_{\tau}}}\cdot$$ =

The values $\tau_{\rm min}$ and $\tau_{\rm max}$ define the range in each ellipsoidal coordinate where the motion is allowed, that is, the coordinate space regions where the condition $p^{2}_{\tau} \geq 0$ is satisfied. The above formulas permit the calculation of the actions $J_{\tau}$ and the frequencies  $\omega_{\tau}$ without integrating the equations of motion. For a Stäckel potential, and given values  E, i₂, i₃ of the integrals of motion, the actions follows from (7), and the frequencies can be computed from (14) and (17).

Thus, in order to compute the frequencies for particular values of H, I₂ and I₃, we have to integrate (17) by numerical means and obtain the roots of $p^{2}_{\tau}=0$, which provides the values of  $\tau_{\rm min}$ and  $\tau_{\rm max}$. This numerical procedure restricts the accuracy with which the frequencies are obtained. We have calculated the integration limits in (17) by means of the zbrent subroutine (Press et al. 1994), setting the tolerance parameter equal to 10^-6. Further, we have performed the numerical integration for (17) by means of the Romberg method -implemented in the midpoint and qromo subroutines (Press et al. 1994) - with an accuracy of the order of 10^-10. Also, for particular values of the integrals H, I₂ and I₃, the functions involved in (17) are not well defined for  $\tau_{\rm min}$ and  $\tau_{\rm max}$. This shortcoming has been avoided by shrinking the integration ranges to ( $\tau_{{\rm min}}+ \delta\tau, \tau_{{\rm max}}- \delta\tau$), ${\rm d}\tau\sim 10^{-4}$ being a suitable value as a compromise between the need to minimize the final errors in the frequencies and the computational time required to the whole numerical procedure.

From the analysis of all the above factors that affect the accuracy with which the frequencies can be obtained, in practice we take as the resonance condition the relation

$$% \vert\bm{m}\cdot \bm{\omega}(\bm {J})\vert\approx 10^{-4},$$ (17)

which determines the different resonant surfaces in the integrals space (H, I₂, I₃).

In order to visualize the resonance structure of the Stäckel model, we fix a value for the energy H, which defines the energy surface in the integrals space. Figure 2 shows the resonance structure for the Perfect Ellipsoid in the I₂I₃-plane, corresponding to H=-0.3, this particular value of the energy being in the range $-0.9074\leq E< 0$corresponding to the Stäckel model, representative of the dynamical behavior of the system. The curves shown in this figure arise from the intersections between the energy surface and several resonant surfaces calculated from (17), in terms of I₂ and I₃, for different resonant vectors m satisfying the condition  |m|=|m₁|+|m₂|+|m₃| < 8. As far a we know, this is the first attempt to compute the resonance structure for this model.

	Figure 2: Resonance structure of the Perfect Ellipsoid onto the I2I3-plane, for a particular value of the energy E = -0.3. The curves shown in this figure are the intersection of the energy surface and several resonant surfaces calculated from (17) for different resonant vectors m satisfying |m|=|m1|+|m2|+|m3| < 8 (see text).
Open with DEXTER

In order to understand how the resonance structure of the Stäckel model determines the orbital structure of the system, in Fig. 3 we present three particular resonances, namely, (1,-3,2),(1,-2,1) and (3,-1,-1), as well as three tori located in different regions on the I₂I₃-plane. Over these tori, orbits will proceed with different dynamical features depending on their location with respect to the resonant curves. Thus, the non-resonant torus depicted in the figure lies in a region of the phase space where the resonant condition does not hold for any m(or where the diophantine condition (11)holds), and any orbit lying on such a 3D torus will cover it densely and uniformly as $t\rightarrow\infty$. As mentioned above, in the 3D torus the motion is ergodic. Meanwhile, motion on a torus located on a resonant curve is such that any trajectory inhabits a submanifold of dimensionality two as a consequence of the resonance condition. Therefore, an orbit on such a torus will densely cover a 2D (elliptic) torus, the resonant torus. Closure in configuration space requires an additional, independent resonance condition. Any torus located on the intersection of two resonances will lead to a periodic orbit since the two independent resonance conditions restrict the motion to one dimension (a 1D torus).

Figures 1 and 2 show the orbital structure in the Stäckel model. The internal curves in Fig. 1, given by I₂=0 and $E=V_{\rm eff}(-\beta)$, play an important role as they bound as well as separate the four regions corresponding to the four families of orbits of the Perfect Ellipsoid. The separatrices are very unstable, even a very small perturbation dramatically modifies the dynamics on them. In the integrable model the smooth separatrices are present because the stable and unstable manifolds of a given object (in this case a 2D hyperbolic torus) exactly match. However, in the surroundings of the separatrices, the resonance structure is rather intricate.

	Figure 3: The three different resonances  (1,-3,2),(1,-2,1) and (3,-1,-1) shown in the figure were extracted from Fig. 2 (see text).
Open with DEXTER

The resonance structure of the integrable Stäckel model allows us to forecast what the dynamical behavior of the system could be if a non-integrable perturbation would be added. As it will be discussed in the next seccion, the interaction of resonances are responsible for the start of chaotic motion when a perturbation is turned on.

   
4 Theoretical discussion of the perturvative effects

Though a deep insight into the dynamics of the perturbed Stäckel system gained by numerical means will be addressed in a forthcoming paper, several of its main features can be inferred from the resonance structure of the integrable model.

We review a theoretical discussion of how a perturbation would modify the dynamics of the integrable Stäckel model, following the general description given by Chirikov (1979) and Cincotta (2002). We consider a system represented by the Hamiltonian

$$% H(\bm{J}, \bm{\theta})= H_{0}(\bm{J}) + \epsilon V(\bm{J},\bm{\theta}),$$ (18)

where

$$% \epsilon V(\bm{J},\bm{\theta})= \epsilon\sum_{\bm{m}} V_{\bm{m}}(\bm {J}) \cos~(\bm{m}\cdot\bm{\theta}).$$ (19)

Here $(\bm {J},~\bm{\theta})$ denote the tridimensional action-angle coordinates for the unpertubed Hamiltonian H₀, m is a tridimensional integer vector, V_m certain real functions and $\epsilon$, the perturbation parameter, is a real number.

Elliptical galaxies might be represented by a smooth triaxial potential H₀, like the Stäckel one, plus a perturbation given by a multipolar expansion, $\epsilon V$. We therefore assume that the perturbation is analytic, which may not be the case if we add a central singularity.

For $\epsilon=0$ we have H=H₀(J) and the motion is stable for any initial condition, since we have the complete set of the three integrals of motion J_i (H₀ being cyclic in  $\bm{\theta}$), i.e. the Hamiltonian is completely integrable. The phase vector evolves linearly with time with a frequency vector $\bm{\omega} (\bm{J})=\partial H_0/\partial\bm{J}$, with ${\rm det}~(\partial\omega_i/\partial J_j)\ne 0$ so that $\bm{\omega}(\bm{J})$ is a one-to-one function. This condition guarantees the nonlinearity of the system.

The presence of the perturbation in general disrupts the integrability of the Hamiltonian H₀ (due to the phase dependence of V), leading to a variation of the unperturbed integrals J_i. The stability of the motion breaks down when a large change in the actions takes place, i.e., when a "gross'' instability sets up.

To describe the motion of a star in the Hamiltonian (18) and (19), we take advantage of a perturbative technique to obtain approximate solutions, assuming that the perturbation parameter is small, $\epsilon\ll 1$, that is, we consider a near-integrable Hamiltonian system. The perturbative approach given by the so-called asymptotic series implies, roughly, that the variation of the unperturbed actions is computed via a power series in the perturbation parameter.

It is well known that the effect of a perturbative Fourier component (as in (19)) is stronger as long as the time variation of its phase, $\dot \psi_{\bm{m}}= \bm{m}\cdot\dot{\bm{\theta}},$ is slow. In the limit case of constant phase we reach the resonance condition for the unperturbed frequencies. If we are far from a resonance (i.e., the initial conditions correspond to values of the actions or the frequencies that satisfy the diophafntine condition (11)), it can be shown that the motion is stable. On taking advantage of the so-called averaging method, we neglect the oscillating part of the perturbation, retaining only its average value.

When we are near to a resonance condition, the asymptotic series technique no longer work due to the appearance of the so-called small denominators in the coefficients of the asymptotic series. These small denominators are the resonant values $\omega_{\bm{m}}=\bm{m}\cdot\bm{\omega}$ that may lead to divergent series. It can be shown that the set  $\omega_{\bm{m}}$ for all integer vectors is, in general, dense everywhere in phase space. Therefore, to find initial conditions "far from a resonance'' is not straight forward.

The geometric features of resonances in the action or integral space are represented in Figs. 2 and 3. Any point on this plane is a torus since the "position vector'' is given by the three values of the actions (or the integrals). In action space, the resonance equation $\bm{m}\cdot\bm{\omega} (\bm{J})=0$leads to another 2-dimensional surface, whose local normal at the resonant point J=J^r is

$$% \bm{n}^r=\left(\partial [\bm{m}\cdot\bm{\omega} (\bm{J})]/ \partial\bm{J}\right)_{\bm{J}^r}.$$ (20)

Further, the conservation of the unperturbed energy provides the 2-dimensional energy surface H₀(J)=E. In what follows we will consider only convex Hamiltonians, a condition that is verified by the Stäckel Hamiltonian. The subspace defined by the intersection of both the resonant and the energy surfaces has dimension 1.

By definition, the frequency vector  $\bm{\omega}$is normal to the energy surface in the action space. The latter condition, together with the resonance equation, shows that the resonant vector m lies in the plane tangent to the energy surface at J=J^rFurthermore, a simple inspection of the equations of motion for only one resonant perturbing term, that is $\epsilon V=\epsilon V_{\bm{m}_g}(\bm {J}) \cos(\bm{m}_g\cdot\bm{\theta})$, shows that  $\dot{\bm{J}}$ is parallel to the resonant vector mg. This picture of the dynamics allows us to conclude that the motion under a single resonant perturbation proceeds on the tangent plane to the energy surface at the point J=J^r along the direction of the resonant vector. As $\epsilon\rightarrow 0$, the resonant perturbation preserves the unperturbed energy.

The above discussion considering only one resonant term shows that the perturbation (19) only depends on a single phase, namely the resonant phase  $\bm{m}_g\cdot\bm{\theta}$. If we perform a canonical transformation in such a way that one of the three new phases, say $\psi_1$, is such that $\psi_1= \bm{m}_g\cdot\bm{\theta}$, then the resulting Hamiltonian will be cyclic in the remainding two phases $\psi_2,~\psi_3$ and the new momenta p₂, p₃will be integrals of the motion.

This canonical transformation may be done through

$$% \psi_i=\mu_{ik}\theta_k;\qquad J_i=J^r_i+p_k\mu_{ki},$$ (21)

where the sum over repeated indexes is understood and $\mu_{ik}$ is a 3 $\times$ 3 matrix. This linear transformation could be thought of as a local change of basis where the new basis vectors are:

$$% \bm{\mu}_{1}=\bm{m}_{g},\qquad\bm{\mu}_{2}={\bm{\omega}(\bm... ...t\bm{\omega}(\bm{J}^r)\vert}}, \qquad \bm{\mu}_{3}\equiv\bm{e}$$ (22)

the being latter normal to both $\bm{\mu}_{2}$ and to n^r. In general e will not be normal to m_g.

According to the discussion given by Chirikov (1979) and Cincotta (2002), the Hamiltonian (18) and (19) in the vicinity of a single resonant perturbation reduces to a simple pendulum model:

$$% H_{r}(p_{1},\psi_{1})\approx H_0({\bm J}^r)+{p_{1}^{2}\over{2M}}+\epsilon V_{\bm{m}_g}(\bm{J}^{r}) \cos\psi_{1},$$ (23)

where H_r stands for the resonant Hamiltonian and M is the "non-linear mass'' given by: $1/M=m_{gi}\left(\partial\omega_i/\partial J_k\right)_{\bm{J}^r}m_{gk}$. Thus p₁ changes with time following the pendulum model for a given value of the integral H_r in the m_g direction. The remaining momenta p_k, k=2, 3are set to zero so that J^r is a point of the trajectory. Therefore, since H_r is independent of time, we again have the full set of three (local) integrals of motion for the problem of a single resonant perturbation: H_r, p₂, p₃.

The variations of the unperturbed actions are determined by the strength of the perturbation and the non-linear mass, in the fashion (see Cincotta 2002 for further details)

$$% (\Delta J)_{r}\equiv\vert\bm{J}-\bm{J}^r\vert\sim \sqrt{\vert V_{\bm{m}_g}(\bm{J}^r)M\vert\epsilon}.$$ (24)

Thus, the resonant surface of the unperturbed system in action space becomes a resonant layer, whose width is given by  $(\Delta J)_{r}$.

Since $\bm{J}-\bm{J}^r\equiv\bm{p}$ and in the new basis $\bm{p}=p_1\bm{\mu}_1+p_2\bm{\mu}_2+p_3\bm{\mu}_3$ we note that p₂ and p₃ being local integrals, the motion proceeds from the resonant point in the direction of  $\bm{\mu}_{1}$, i.e. across the resonance.

Cincotta (2002), following Chirikov (1979), discusses how different terms in the pertrubation modify the dynamics of the integrable system H₀. For clarity we recast (18) and (19) as

 

$$% H(\bm{J}, \bm{\theta}) H_{0}(\bm{J}) + \epsilon V_{\bm{m}_g}(\bm {J}) \cos~ (\bm{m}_g\cdot\bm{\theta})$$ =

$$+ \epsilon V_{\bm{m}_l}(\bm {J}) \cos~ (\bm{m}_l\cdot\bm{\theta})... ...ne \bm{m}_g, \bm{m}_l}\!\! V_{\bm{m}}(\bm {J})\cos~(\bm{m}\cdot\bm{\theta}).\ ~$$ (25)

If we take into account all terms in (25), all momenta change with time and then variations in p₂ and p₃ could occur, since now the Hamiltonian depends on  $\bm{\theta}_2$ and  $\bm{\theta}_3$. These variations may proceed along the  $\bm{\mu}_{2}$ and $\bm{\mu}_{3}$ directions, and recalling that  $\bm{\mu}_{2}$ is normal to the unpertrubed energy surface, the variations in p₂ take into account small changes in the unperturbed energy when the pertubation is switched on. These variations are bounded and of order $\epsilon$. However, $\bm{\mu}_{3}$ is tangent to the energy surface and by definition it is directed along the resonance, so one could expect variations in p₃ that, at first sight, will be unbounded (see next section).

Let us now discuss the effects of the different terms in (25). As shown above, on taking V_m=0 for all $\bm{m}\ne\bm{m}_g$and picking up initial conditions such that $\bm{m}_g\cdot\bm{\omega}(\bm{J}^r)=0$, the torus structure in the vicinity of J^r is preserved because of the existence of the three local integrals H_r, p₂ and p₃, $H_r(p_1,\psi_1)$ being the pendulum Hamiltonian. Depending on the value of H_r, the momenta p₁ will oscillate or rotate around J^r. The resonant layer corresponds to the oscillation domain of the pendulum energy. As it is well-known, oscilations and rotations in a pendulum are separated by a smooth separatix  $p_1^s(\psi_1^s)$.

If we switch on V_ml, which we assume to be the dominant non-resonant perturbing term, its main effect on the dynamics of the system in the vicinity of J^r would be to split the stable and unsable manifolds associated with the unstable points of the pendulum leading to a splitting of separatices and the smooth curve  $p_1^s(\psi_1^s)$becomes a stochastic layer of finite width. Thus, under the effect of this perturvative term, oscillations and rotations are actually separated by a layer of chaotic motion; p₁(t) shows a bounded unstable, chaotic behavior confined to this thin layer. Figure 4 illustrates this fact for two given resonances in the Stäckel model in the unperturbed integral space  (H₀,I₂,I₃) for a fixed value of H₀. There, the local basis $\{\bm{\mu}'_{1}, \bm{\mu}'_{2}, \bm{\mu}'_{3}\}$, at a given resonant torus I^r, has been included, the  $\bm{\mu}'_{i}$ the  $\bm{\mu}_{i}$ at J^r being mapped to the integral space.

	Figure 4: The upper panel shows the resonances (1,-3,2) and (3,-1,-1) of Fig. 3 (reproduced in left panel), with their concomitant resonant and stochastic layers.
Open with DEXTER

If we turn on the next relevant term in (25), say V_md, it is possible to show that this third term might cause unbounded chaotic variations in p₃while p₁ lies in the stochastic layer; there could be motion in the $\bm{\mu}_{3}$ direction, that is, along the stochastic layer of the resonance.

Figure 5 illustrates in a schematical way how the whole resonant structure of the Stäckel model (given in Fig. 2) could be distorted by the effect of a generic perturbation. All resonaces would become layer, whose widths  $(\Delta J)_{r}$strongly depend on J^r, m, V_m but all of them have order $\epsilon$.

	Figure 5: Schematic representation of the effect of a perturbation on the resonant structure of the Stäckel model (given in Fig. 2). All the resonances become layers of width  $(\Delta J)_{r}$.
Open with DEXTER

Geometrically, when the system is exactly at some resonance, we have a 2D elliptical torus, while when the system is on the border of the resonance we get a 2D hyperbolic torus. Thus the center of a resonance layer over the integral space form a chain of elliptical tori while the border of the resonace layer corresponds to a chain of hyperbolic tori. This is schematically shown in Fig. 6 for the same resonances of the Stäckel potential depicted in Fig. 5.

	Figure 6: Geometric structure of the perturbed resonances shown in Fig. 5. The central curve, that corresponds to the exact resonance is a chain of 2D elliptical tori, while the borders of the resonance correspond to a chain of 2D hyperbolic tori.
Open with DEXTER

All the orbits trapped in resonances constitute subfamilies of the Stäckel orbital families. As we discussed at the beginning of this section, if the system is far from any resonance, the motion is stable. The presence of the perturbation only produces small changes of order $\epsilon$ in the unperturbed actions. Thus a non-resonant value of the unperturbed actions or integrals leads to a 3D non-resonant torus that support quasiperiodic orbits. Accordingly to Eqs. (7) and (8), the three local integrals are $H=H_0+\epsilon V,~ I_2'=I_{2}+\epsilon\delta I_2,~ I_3'=I_{3}+\epsilon\delta I_3$, $\delta I_i$ being bounded constants. Therefore, in the neighbourhood of the values of the unperturbed integrals ( I₂,I₃) the orbital structure is preserved.

On the other hand, orbits trapped in a single resonance at J^r also respect three local integrals that in action space are the pendulum Hamiltonian H_r and two independent linear combinations of the unperturbed actions at the resonant value, for instance, K₁=a₁J₁^r+a₂J₂^r,  K₃=a₃J₁^r+a₄J₃^r, where the coeficients a_i depend only on the $\mu_{ki}$. For H_r=0 we are exactly on a 2D elliptical torus that supports a stable resonant orbit that is the same as in the unperturbed system (since $p_1\equiv 0$and J=J^r). For positive values of H_r new 3D tori appear that support nearly resonant orbits that respect the above mentioned local integrals. The nature of these new orbits are different to that of the unperturbed system at J^r. When H_r takes the value that corresponds to the pendulum separatrix, we arrive at to the 2D hyperbolic tori that support unstable orbits.

   
5 Discussion on the origin of chaos and diffusion

The dynamical description becomes highly intricate when we consider the interaction among all the perturvative terms. Indeed, as mentioned above, the effect of a perturbation on resonance is to form a thin stochastic layer at the border of the resonance layer. Since the resonace layer has a width $\sim$ $\! \sqrt{\epsilon}$, as long as the perturbation increases its strength, the resonaces become wider and the overlap of the stochastic layers of different resonances coud take place when $\epsilon$ reaches some critical value  $\epsilon_c$. When this occurs, the resonances overlap. The dynamics in this case become unstable due to the intersections of the stable and unstable manifolds of different resonances. The resulting motion is completely stochastic (see Cincotta 2002) whose main characteristic is its local instability. Figure 5 provides a schematic illustration of the effect of a perturbation on the resonant structure. From this figure. we can clearly see that if we increase the perturbation, a massive overlap of resonaces will occur. Due to the dimensionality of the system (3D), resonance intesections take place (which is not the case in 2D systems). Around the point of intersection, the stochastic layers of different resonaces ovelap, as Fig. 5 shows schematically.

For the time being, an analytical description of the dynamics at the resonance intersection is still lacking, because if we retain two resonant terms in the expansion (25), for instance V_mg and  V_ml, the Hamilonian depends then on two resonant phases  $(\bm{m}_g\cdot\theta)$ and  $(\bm{m}_l\cdot\theta)$, and as can be shown, this system is not integrable. A picture of the dynamics at a resonance crossing is given in Cincotta et al. (2003), where a zoom of the dynamical behaviour at the region close to the intersection of two resonances with their corresponding stochastic layers is presented. For this case, the intersection of these two resonances generates a stable periodic orbit and a central region of regular or mild chaotic motion bounded by a thin layer of unstable motion that presents at the same time several stable domains separated by extremely thin chaotic filaments. However, it is possible to infer that at the intersection of the unperturbed resonances a periodic orbit appears whose stability cannot be determined a priori and that the zone of stochastic motion should be bounded (see Chirikov 1979).

The kown mechanisms that lead to a transition from regularity to chaos are overlap of resonances (including resonace crossings) and Arnold diffusion-like processes (see for instance, Chirikov 1979; Cincotta 2002; Giordano & Cincotta 2004).

As we have already stated, chaos means variation of the unperturbed integrals. This is usually called chaotic diffusion in the literature. Unfortunatelly no theory exists that describes global diffusion in phase space. In other words, it is not possible to estimate either its rate or its direction or route. Although one can obtain accurate values of the Lyapunov exponents, the KS entropy, the MEGNO (see Cincotta et al. 2003; Cincotta & Simó 2000) or any other indicator of the stability of the motion, they provide only local values of the diffusion rate. A given orbit in a chaotic component of the phase space coud have two positive and large values for the Lyapunov exponents which does not mean that the unperturbed integrals would change much. This is a natural consequence of the structure of the phase space of almost all dynamcal systems like galaxies.

Thus what is actually relevant is the extent of the domain and the time-scale over which the diffusion may occur. In a recent work Giordano & Cincotta (2004) showed numerically that in models similar to that of an elliptical galaxy, the time-scale over which diffusion becomes relevant is several orders of magnitude of the Hubble time. On the other hand, in models corresponding to planetary or asteroid dynamics, diffusion may occur over physical time-scales.

One can make predictions about the global dynamical behaviour of a perturbed system by studying its resonance structure. For example, in Fig. 2 we observe a rather intrincate resonance structure in the box orbits domain, with several crosses of resonance as well as the fact that they are very close to each other. Thus, we expect that this region would be dominated by chaotic dynamics when the perturbation is actived, for instance by adding a mass concentration in the center of the system, the box domain in the I₂I₃-plane will be the first to be dominated by chaos. This is also true for any other kind of integrable model and perturbation, whatever its strengh at the center may be (see, for instance Papaphilippou & Laskar 1998).

We could come to a similar conclusion for those regions close to the separatrices which separate different families of orbits, I₂=0 and $E=V_{\rm eff}(-\beta)$. Several resonance intersections along these separatices can be observed, then, when a non-integrable perturbation is added, the dynamics in this region would become rather intrincate due to the fact that the stable and unstable manifolds of different 2D hyperbolic tori will tangle in a very complicate fashion (see, for instance, Cincotta et al. 2003). This explains the instabilities found by de Zeeuw (1985) along the I₃-axis; in fact, in the Stäckel model all 2D tori with I₂=0 are unstable (hyperbolic), so that the dynamics in this region becomes rather chaotic under the effect of a non-integrable perturbation.

An important fact to be stated is that when chaos sets up, the unperturbed global integrals (or actions) have a discontinuous dependence on phase space variables. Indeed, close to the resonant torus, despite the existence of three local integrals, the unperturbed orbital structure is not preserved and the topology of the phase space changes. Moreover, on the stochastic layer at least one integral does not exist.

Thus, in the Stäckel model, for $\epsilon=0$ we get Fig. 2, while for $\epsilon\ll\epsilon_c$ the dynamic picture will be similar to that given in Fig. 5. Resonances do not present a significant overlap and chaotic motion will be confined to the narrow stochastic layers and around the ressonance crossings. Thus we cannot expect large variations of the unperturbed integrals and therefore chaotic diffusion should be irrelevant. However, theoretically, it might be possible that chaotic diffusion takes place even for very small values of the perturbation ( $\epsilon\ll\epsilon_c$). Since for very small $\epsilon$ values the resonance structure is still present, one could expect that the perturbative terms drive diffusion along the stochastic layers of a given resonance, i.e. in the $\bm{\mu}'_{3}$ direction.

Indeed, Arnol'd (1964) proved for a very specific Hamiltonian, that there exists a mechanism that could connect two distant tori in the chain of 2D hyperbolic tori located on the border of the resonance. This means that variations of the unperturbed actions along the resonance occur. Since the actual motion around any hyperbolic tori is stochastic, chaotic diffusion along the stochastic layer of a given resonance could take place. The mechanism that allows motion in the $\bm{\mu}_{3}$ direction is the so-called Arnol'd Mechanism (see Giorgilli 1990). However, the existence of a transition chain of hyperbolic tori does not mean that stochastic motion would spread over the whole web of resonances. Moreover, in general, it is not possible to find this transition chain in most non-linear Hailtonian systems, even in the case of the simpler ones. Nevertheless, it is often assumed that this is the case and that chaotic motion will take place over the whole resonance web. This implies that a global instability exists and it is usually called Arnol'd diffusion. This conjecture supports, for instance, concepts like mixing or ergodicity, among others. Let us recall that the global instability properties of near integrable Hamiltonian systems are far from being well-understood after the pioneering work of Arnol'd (see Lochak 1999).

For $\epsilon\gtrsim\epsilon_c$, the overlap of resonances begins to operate destroying the torus structure of phase space and the resonance web disappears, and chaotic diffusion begins. In this case, we expect that diffusion could operate through overlap of resonances, but the extent of the diffusive zone will be confined, at most, to the overlapping domain. There is no guarantee that the whole chaotic component of phase space is connected, even in the case of massive overlap (Giordano & Cincotta 2004). But also in this case, it is usually stated that the chaotic component of phase space is fully connected through Arnol'd diffusion. Arguments based on these yet unproved statements have given rise to many of the current ideas taken for granted on dynamics in elliptical galaxies, as can be verified for instance in the works of Merritt & Valluri (1996), Merritt (1999) and Udy & Penniger (1988) among others.

Therefore, in a theoretical approach, while it is posible to get much information about the dynamical behaviour of the system as long as the perturbation increases its strength, nothing can be concluded about the diffusive processes in phase space of near-integrable systems. Numerical integrations become the only way to investigate if chaotic diffusion could play a significant role in particular models of galaxies. This will be addressed in a forthcoming paper using the Stäckel potential discussed here.

We clearly discriminate the dynamical behaviour of the system depending on the values of the integrals. Close to strong non-resonant tori, the orbital structure of the unperturbed Stäckel potential is preserved and the local integrals are corrections of order $\epsilon$of the unperturbed global integrals. On the other hand, when the system is close to a resonant tori, the unperturbed orbital structure is not preserved, new local integrals appear and the topology of phase space changes. These new local integrals are the pendulum Hamiltonan and linear combinations of the unperturbed actions at the resonant point.

Therefore, in such a system, it is not possible to assume that the distribution function of the galaxy, in the whole regular component, has the form  f(I₂, I₃). This could be true only for strong non-resonant tori, but since resonances are dense in phase space, the distribution function locally exists and has the form  $f_n(H, I_{2}, I_{3})+\epsilon g(H, I_{2}, I_{3})$ in neighbourhood of non-resonant tori that support quasiperiodic orbits and  f_r(H_r, K₂, K₃) in vicinity of resonant tori.

On the other hand, nothing can be said about the dependence of f in the chaotic region. In a similar fashion, since there is no theoretical support to argue that the whole chaotic region is fully connected, one should define a local distribution function only in those regions of phase space that could be visited by a single orbit. It is very likely that in some regions two local quasi-invariants exist (for example in the stochastic layers of the resonances) and in some other regions only the energy remains constant. Clearly, discontinuous dependence of f on the integrals, and consequently on phase space variables, is expected.

On the other hand, much progress can be made in the understanding of diffusive processes in phase space of near-integrable Hamiltonian systems, including the dynamical behaviour around resonance crossings. Indeed, as all the evidence shows, galaxies should exhibit a divided phase space, with both regular and chaotic components, which is the phase space structure of near-integrable Hamiltonians for low-to-moderate perturbations.

Taking into account these arguments that are supported by theoretical considerations, we suggest that the direction to take in constructing equilibrium models of elliptical galaxies.

The selected representation of an elliptical galaxy does not encompass the case of a central singularity. In such a case, one deals with a non-analytic perturbation which admits no Fourier expansion as (18) and the perturbative approach is no longer suitable. Nevertheless, it is clear that such a peturbation $\sim$1/r mainly affects the dynamics in the box domain due to the fact that these orbits pass arbitrarily close to the center: box orbits are strongly scattered producing in general an almost spherical configuration near the center (see however Poon & Merritt 2002). Therefore, the perturbative approach addressed here should be reformulated so that a singularity at the origin can be taken into account. This is not a simple task. Although there is observational evidence of central mass concentrations in elliptical galaxies that could be black holes (see for instance Ferrarese & Ford 2005), the question of whether they actually are supermassive black holes is still an open matter.

Due to the fact that the described perturbative approach succeeds in studying actual physical problems such as planetary and asteroidal systems or dynamics in particle accelerators as well as the motion of charged particles in magnetic bottles, we argue that the theoretical formulation discussed here, which the only one developed so far, is a fairly useful tool when dealing with galactic systems.

Although one may argue that the action-angular variables do not constitute an adequate set to yield a description of galatic systems, there exists numerical evidence of galatic models exhibiting large regions of phase space corresponding to regular motion. Therefore, the Hamiltonian of a star moving in a galactic potential can always be writen in the fashion $H(\bm{J}) = H_{0}(\bm{J}) + \epsilon V(\bm{J},\bm{\theta})$, the key point being what should be considered as $V(\bm{J},\bm{\theta})$. Since all these mathematical theories have served well to provide a thorough description of chaotic systems, they may be of use to yield a good first order approximation of the motion of a star in a smooth stationary potential. Moreover, despite the fact that the actual dynamics of a galactic system comes from numerical simulations, theory that permits us to elucidate the matter, at least in an heuristic way is still useful.

Acknowledgements

This work was partially supported by grants of the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) and Fundación Antorchas. The authors aknowledge the valuable comments and suggestions of the referee.

References

• Arnol'd, V. I. 1964, Sov. Math.-Dokl., 5, 581 (In the text)
• Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium (New Haven: Yale University Press) (In the text)
• Chirikov, B. 1979, Phys. Rep., 52, 63 (In the text)
• Cincotta, P. M. 2002, New AR, 46, 13 [NASA ADS] (In the text)
• Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [EDP Sciences] [NASA ADS] (In the text)
• Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Phys. D, 182, 151 [CrossRef] (In the text)
• Crane, R., Stiavelly, M., King, I. R., et al. 1993, AJ, 106, 1371 [NASA ADS] [CrossRef] (In the text)
• de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] (In the text)
• de Zeeuw, T., & Lynden-Bell, D. 1985, MNRAS, 215, 713 [NASA ADS] (In the text)
• Ferrarese, L., van den Bosch, F. C., Ford, H. C., Jaffe, W., & O' Connell, R. W. 1994, AJ, 108, 1598 [NASA ADS] [CrossRef] (In the text)
• Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9 [NASA ADS] [CrossRef] (In the text)
• Ferrarese, L., & Ford, H. 2005, Space Sci. Rev., Vol. 116, 523 (In the text)
• Ford, H. C., et al. 1998, The central Regions of the Galaxy and Galaxies, ed. Y. Sofue (Dordrecht: Kluwer), IAU Symp., 184, 377 (In the text)
• Gerhard, O. E., & Binney, J. J. 1985, MNRAS, 216, 467 [NASA ADS] (In the text)
• Giordano, C. M., & Cincotta, P. M. 2004, A&A, 423, 745 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Giorgili, A. 1990, in Les Métodes Modernes de la Méchanique Céleste, ed. C. Froeschlé, & C. Benest, 249 (In the text)
• Kormendy, J., & Richstone, D. O. 1995, ARA&A, 33, 581 [NASA ADS] (In the text)
• Kuzmin, G. G. 1956, AZh, 33, 27 (In the text)
• Kuzmin, G. G. 1973, in The Dynamics of Galaxies and Star Clusters, ed. T. B. Omarov, 71 (In the text)
• Laskar, J. 1993, Phys. D, 67, 257 [CrossRef] (In the text)
• Lochak, P. 1999, in Hamiltonian Systems with Three or More Degress of Freedom, ed. C. Simó, NATO ASI (Dordrecht: Kluwer), 168 (In the text)
• Merritt, D. 1999, PASP, 111, 129 [NASA ADS] [CrossRef] (In the text)
• Merritt, D., & Ferrarese, L. 2001, MNRAS, 320, L30 [NASA ADS] [CrossRef] (In the text)
• Merritt, D., & Fridman, T. 1996, ApJ. 460, 136 (In the text)
• Merritt, D., & Quinlan, G. 1998, ApJ, 498, 625 [NASA ADS] [CrossRef] (In the text)
• Merritt, D., & Valluri, M. 1996, ApJ, 471, 82 [NASA ADS] [CrossRef] (In the text)
• Papaphilippou, Y., & Laskar, J. 1998, A&A, 329, 451 [NASA ADS] (In the text)
• Poon, M., & Merritt, D. 2002, ApJ, 586, 89 [CrossRef] (In the text)
• Press, W., et al. 1994, Numerical Recipes in Fortran (Cambridge University Press) (In the text)
• Peebles, P. J. E. 1972, Gen. Relativ. Gravitation, 3, 63
• Schwarzschild, K. 1981, The Structure and Evolution of Normal Galaxies (Cambridge University Press) (In the text)
• Schwarzschild, M. 1993, ApJ, 409, 563 [NASA ADS] [CrossRef] (In the text)
• Udry, S., & Pfenniger, D. 1988, A&A, 198, 135 [NASA ADS] (In the text)
• Valluri, M., & Merritt, D. 1998, ApJ, 506, 686 [NASA ADS] [CrossRef] (In the text)
• Wachlin, F., & Ferraz-Mello, S. 1998, MNRAS, 298, 22 [NASA ADS] [CrossRef] (In the text)