A&A 423, 745-753 (2004)
DOI: 10.1051/0004-6361:20040153

Chaotic diffusion of orbits in systems with divided phase space

C. M. Giordano - P. M. Cincotta

Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque, 1900 La Plata, Argentina

Received 12 August 2003 / Accepted 7 April 2004

Abstract
In this paper we discuss the relevance of diffusive processes in multidimensional Hamiltonian systems. By means of a rather simple model, we present evidence that for moderate-to-strong chaotic systems the stochastic motion remains confined to disjoint domains on the energy surface, at least for mild motion times. We show that only for extremely large timescales and for rather large perturbations, does the chaotic component appear almost fully connected through the relics of the resonance structure. The discussion whether diffusion over the energy surface could actually occur in asteroidal or galaxy dynamics is also included.

Key words: diffusion - stellar dynamics - celestial mechanics

  
1 Introduction

In order to build up an equilibrium model for a galaxy that is assumed to be represented by a smooth gravitational field it is necessary to know beforehand its global dynamics. It seems likely that any realistic model should exhibit a divided phase space, that is, the motion would take place either in a stable, regular component or in one or more unstable, chaotic components. In 3D systems displaying such a dynamics, the existence of three, in general local, invariants, allows the presence of invariant tori where regular, quasi-periodic or resonant motion takes place. The disruption of these local invariants, mainly due to resonance interactions, leads to the appearance of a chaotic component.

However, the first attempts to investigate these matters in astrophysical literature assumed that for dynamical systems with more than 2 degrees of freedom the chaotic component is fully connected. In fairly recent studies, such as Merritt & Valluri (1996) and Merritt & Fridman (1996), the authors study whether this full connection of the chaotic component may occur in realistic physical times. If such were the case, it would imply that the orbit of a star within the chaotic domain would explore the whole region, which, in general, would comprise a large fraction of the energy surface. As a consequence the distribution function on the chaotic component would depend only on the energy (see Merritt 1999 for a thorough discussion).

Here, by means of a very simple model, we present numerical evidence that for moderate-to-strong chaotic systems diffusion does not occur over the whole chaotic component and only when the latter fills almost all the energy surface may diffusion become significant but in extremely long times. As far as we know, this is one of the first efforts showing diffusive processes in phase space of multidimensional systems. With this aim we have chosen a rather simple Hamiltonian system in order to have a complete knowledge of the local and global dynamics and the transition to chaos and particularly to be able to study in detail the motion in the resonance intersections.

In a forthcoming paper we will present a detailed investigation of these matters but we will now deal with a fairly realistic model, namely, a perfect ellipsoid perturbed by Denhen's law.

  
2 The resonance structure for a 3D toy model

As is well known, the knowledge of the resonance structure is of actual relevance when dealing with near-integrable Hamiltonian systems.

In the present effort, we will be concerned with the perturbed uncoupled quartic oscillator

$$\tilde{H}(\bm{p},\bm{q})={\bm{p}^2\over 2}+{1\over 4}(x^4+y^4+z^4)+\epsilon x^2(y+z).$$ (1)

The full dynamics of this 3D system has been investigated by Cincotta & Giordano (2000) and Cincotta et al. (2003) by means of the Mean Exponential Growth factor of Nearby Orbits (MEGNO) (for a detailed formulation of this technique see Cincotta et al. 2003 and Cincotta & Simó 2000; a brief description is given in Appendix A). The advantage of this toy model is that it can be easily written in terms of action variables and that the coordinates admit of a simple Fourier expansion. Indeed, in terms of the unperturbed action-angle variables, $(I_1, I_2, I_3; \theta_1, \theta_2, \theta_3)$, the Hamiltonian (1) can be recast as

$$H(\bm{I},\bm{\theta})=H_0(\bm{I})+\epsilon V(\bm{I},\bm{\theta}),$$ (2)

where H₀ is given by

$$H_0(\bm{I})=A\;\big(I_1^{4/3}+I_2^{4/3}+I_3^{4/3}\big),$$ (3)

with $A=(3\beta/2\sqrt{2})^{4/3}$, where $\beta=\pi/2K(1/\sqrt{2})$, K(k) denoting the complete elliptic integral. The perturbation, which can be assumed to be small as long as  $\epsilon\ll 1$, admits of the Fourier expansion

 

$$V\!(\bm{I},\bm{\theta}) \hat{V}_{12} \sum_{n,m,k=1}^{\infty} {\alpha_{nmk}\Big( \cos\big(2(n + m - 1)\theta_1\pm(2k - 1)\theta_2\big)}$$ =

$$+ \cos\big(2(n - m)\theta_1\pm(2k - 1)\theta_2\big)\Big)$$

$$+\hat{V}_{13} \sum_{n,m,k=1}^{\infty} {\alpha_{nmk}\Big( \cos\big(2(n + m - 1)\theta_1\pm(2k - 1)\theta_3\big)}$$

$$+ \cos\big(2(n - m)\theta_1\pm(2k - 1)\theta_3\big)\Big),$$ (4)

the function $\hat{V}_{ij}$ and the coefficients  $\alpha_{nmk}$ being

$$\begin{eqnarray*}\hat{V}_{ij}&=&{3\beta\over 4} I_{\rm i}^{2/3}I_j^{1/3}, \;\;\;... ...;\quad {\alpha_{s+1}\over\alpha_s}\approx {1\over 23},\nonumber \end{eqnarray*}$$

and can be split in two, namely V_xy and V_xz, on introducing the integer vectors  l=(l₁,l₂,0) and  k=(k₁,0,k₃), and the new coefficients  $\hat{\alpha}_{l_1l_2}$ and  $\hat{\alpha}_{k_1k_3}$,

$$\begin{eqnarray*}V_{xy}(I_1,I_2;\theta_1,\theta_2) &=& \hat{V}_{12}\sum_{l_1,l_2... ..._3}\hat{\alpha}_{k_1k_3} \cos(k_1\theta_1+k_3\theta_3).\nonumber \end{eqnarray*}$$

The concomitant unperturbed frequency vector is given by

$$\bm{\omega}(\bm{I})={\partial H_0\over\partial\bm{I}}={4\over 3}A\big(I_1^{1/3},I_2^{1/3},I_3^{1/3}\big).$$ (5)

For the unperturbed Hamiltonian H₀, the resonance structure on the energy surface can be easily visualized by introducing a change of coordinates such that the unperturbed energies in each degree of freedom, h₁, h₂, h₃, become the new action-like variables. Therefore we obtain

$$\bm{\omega}(h_1,h_2,h_3) = \sqrt{2}\beta \big(h_1^{1/4},h_2^{1/4},h_3^{1/4}\big),$$

and, in terms of (h₁, h₂), the resonance condition, $\bm{m}\cdot\bm{\omega}=0,~ \bm{m}\in{\mathchoice {\hbox{$\sf\textstyle Z\kern-... ...kern-0.3em Z$ }} {\hbox{$\sf\scriptscriptstyle Z\kern-0.2em Z$ }}}^3/\{\bm{0}\}$, for H₀=h can be recast in the form

 

$$(m_1^4+m_3^4)\xi^4 + 4m_1^3m_2\xi^3\eta + 6m_1^2m_2^2\xi^2\eta^2 + 4m_1m_2^3\xi\eta^3 + (m_2^4+m_3^4)\eta^4 - m_3^4=0,$$ (6)

where $\xi=(h_1/h)^{1/4},~ \eta=(h_2/h)^{1/4}$.

Note that for those harmonics in which one of the $m_{\rm i}$ is zero, the resonant polynomial (6) can be easily solved to yield

   

$$\hat{h}_2 {m_1^4\over m_2^4}\hat{h}_1,\quad m_1m_2<0, \quad m_3=0,$$ = (7a)

$$\hat{h}_2 1-\Big({m_1^4+m_3^4\over m_3^4}\Big)\hat{h}_1, \quad m_1m_3<0, \quad m_2=0,$$ = (7b)

$$\hat{h}_2 \Big({m_3^4\over m_2^4+m_3^4}\Big)(1-\hat{h}_1),\quad m_2m_3<0, \quad m_1=0,$$ = (7c)

with $\hat{h}_{\rm i}=h_{\rm i}/h$, showing that those resonances associated to resonant vectors with at least one null $m_{\rm i}$ appear as straight lines on the energy surface $\hat{h}_1+\hat{h}_2+\hat{h}_3=1$.

The width of any of those resonances can be computed by means of a simple pendulum approximation, which indeed is a suitable description whenever we assume each resonance to be isolated from the rest.

Let us notice that the map ${\vec I} \mapsto {\vec h}$ transforms the unperturbed energy surface into a plane, so that we can perform a second (global) change of coordinates, $(h_1,h_2,h_3)\mapsto (e_1,e_2,e_3)$, in such a way that the e₃-axis is normal to the energy plane,

   

$${1\over\sqrt{6}}(h_1-2h_2+h_3),$$                           e₁ = (8a)

$${1\over\sqrt{2}}(h_1-h_3),$$ e₂ = (8b)

$${1\over\sqrt{3}}(h_1+h_2+h_3),$$ e₃ = (8c)

with

 

$$-\sqrt{2\over 3}\le {e_1\over h}\le {1\over\sqrt{6}},\quad -{1\ov... ...}\le {e_2\over h}\le {1\over\sqrt{2}},\quad {e_3\over h}= {1\over\sqrt{3}}\cdot$$ (9)

The resonance structure of the unperturbed system on the energy surface and the theoretical widths of the principal resonances appearing in the perturbation are presented in Figs. 1a and 1b, respectively.

In Fig. 1b, the resonance widths were measured in terms of  $\bm{\Delta h}$ instead of  $\bm{\Delta I}$ and then transformed to the (e₁,e₂)-plane by means of Eqs. (8); for the detailed computation of the resonance widths we refer to Cincotta et al. (2003).

On applying the MEGNO (a somewhat detailed description of this tool can be found in Appendix A) to the study of the global dynamics of Eq. (1), one obtains the resonance structure presented in Fig. 1c. There, the details of the phase space structure at a low-to-moderate value of the perturbation ( $\epsilon =5\times 10^{-3}$) are displayed, depicting the obtained values for the MEGNO in a contour-like plot where resonances can be clearly distinguished.

Figure 1: a) Resonances of the unperturbed Hamiltonian (3) for 1<|m1|+|m2|+|m3|<9 and energy h=1yielding the theoretical Arnold web on the energy surface. b) Strongest resonances and their theoretical widths. Arrows within resonances (2,-3,0) and (2,0,-3) indicate the direction in which $\bm {\Delta h}_{\bm {m}}^r$ oscillates about the corresponding resonant value. (Both Figs.  a) and  b) have been taken from Cincotta et al. 2003.) c) Actual dynamics of the system revealed by the MEGNO: $\overline {Y}(t_{\rm f})$-levels on the energy surface for $\epsilon =5\times 10^{-3}$, $h=1/(4\beta ^4)\approx 0.485$ (see text). The contour plot corresponds to $\overline {Y}$ binned in the intervals [1.99, 1.995), [2, 2.015),  [2.015, 20),  [20,  160),  [160,  215) , being $t_{\rm f}=3500T$.

Open with DEXTER

At this stage, a brief reference to MEGNO's behavior is required in order to grasp the dynamical information comprised in such a figure. Let us then recall that in the case of regular motion, the MEGNO ( $\overline {Y}$) tends asymptotically to a fixed value independent of the orbit, namely $\overline Y \to 2$. Small departures from this value indicate the proximity of some periodic orbit, where $\overline Y \la 2$ and  $\overline Y \ga 2$ for stable and unstable periodic orbits, respectively. On the other hand, for irregular, stochastic motion  $\overline {Y}$ grows linearly with time as  $t \to \infty$, at a rate equal to $\sigma/2$, $\sigma$ denoting the largest Lyapunov Characteristic Number (LCN) of the orbit.

Let us devote the following paragraphs to describing the procedure giving rise to Fig. 1c, and Figs. 2a and 2b which correspond to higher values of the perturbation.

	Figure 2: $\overline {Y}(t_{\rm f})$-levels on the energy surface. The contour plot in gray scale (from white to dark-gray) corresponds to  $\overline {Y}$ binned in the intervals, a)  [0, 2.01),  [2.01, 15),  [15, 50),  [50, 160),  [160, 412); b) [0, 2.1),  [2.1, 30),  [30, 60),  [60, 130),  [130,310)  [310, 512).
Open with DEXTER

Though at first sight the Hamiltonian (1) depends on two parameters - the energy h and $\epsilon$ -, after rescaling the variables in the fashion  $x_{\rm i}\to\epsilon x_{\rm i},~ p_{\rm i}\to\epsilon^2p_{\rm i}$ and $t\to t/\epsilon$, we have  $\tilde{H}\to\epsilon^4\tilde{H}\equiv\bar{H}$, which is independent of $\epsilon$, and therefore the scaled energy, $\bar{h}=\epsilon^4 h$, turns out to be the only free parameter. We then allow $\epsilon$ to vary and fix the energy at the value $h=1/(4\beta ^4)\approx 0.485$. This adopted value for the energy leads to a period $T=2\pi$ for the y, z-axial periodic orbits, which remain always stable despite the strength of the perturbation.

For each adopted value of $\epsilon$, we take values of h₁ and h₂ with $0\le h_1,h_2\le h,~ h_3=h-h_1-h_2$, where h₁ and h₂ are of the form  $jh/1000;~j=0,\ldots,1000.$ This leads to 501, 501 initial conditions for which we take (x,y,z)=(0,0,0). We integrate the equations of motion together with their first variationals over a total motion time $t_{\rm f}=3500T$. For the tangent vector, we adopt the initial values $\delta_x=\delta_y=\delta_z=0$ and  $\delta_{p_{\rm i}}$ chosen at random in the interval (-1,1) and then normalized to 1. For each orbit we compute both  $\overline {Y}(t_{\rm f})$ and the rate at which the MEGNO grows with time. When performing the least squares fit on  $\overline{Y}(t)$ to obtain the LCN, only the last $80\%$ of the time interval is considered in order to avoid the initial transient. The actual energies  h₁,h₂, h₃, are scaled to the interval [0,1] (by division through h) and then transformed to the energy plane  (e₁, e₂) by means of Eqs. (8).

In Fig. 1c, corresponding to $\epsilon =5\times 10^{-3}$, we display the values of  $\overline {Y}(t_{\rm f})$ binned in five intervals, two of them being very narrow and close to 2 (see figure caption for details). We have skipped all points with $1.995\le\overline{Y}(t_{\rm f})<2$(corresponding to quasi-periodic motion) so that the resonances can be discerned more clearly. The gray-scale has been simulated using different point sizes for different  $\overline {Y}(t_{\rm f})$ intervals, which have been selected such as to highlight some details of the dynamics at the chosen perturbation value.

In this picture most of the main resonances can be clearly distinguished as light gray channels surrounded by dark boundaries. Thus, for instance, four main resonances are seen to intersect at the origin, the three lines corresponding to the  (1,-1,0),(1,0,-1),(0,1,-1) resonances and the curve associated to the (-2,1,1) resonance.

The actual width of the resonances as well as the narrow stochastic layers at their edges can be clearly visualized. The center of any resonance "channel" corresponds to a sequence of 2D elliptic tori while its borders ("stochastic layer") correspond to a sequence of 2D hyperbolic tori. We observe a strip of chaotic motion close to h₁=0. The presence of this region is easily understood from Fig. 1a as the overlap of the strongest resonances, e.g. (2,-1,0) and (2,0,-1), as well as of many others (not shown in Fig. 1a).

It is quite interesting to see how intricate the dynamics at the intersection of resonances can be. To illustrate such a feature, let us note the complexity of the picture reproduced in Fig. 5c, taken from Cincotta et al. (2003), where we present a zoom around the intersection of resonances at the origin in Fig. 1c. That contour plot was obtained with a higher resolution in h₁ and h₂ and for a total motion time  $t_{\rm f}=350T$. There the MEGNO reveals the existence of several stability zones, which should be responsible for restraining the spread of chaotic motion, acting as barriers to diffusion. They are the well-known sticky tori surrounding the periodic orbit located at the center of the resonance. This plot is also very illustrative to see how the manifolds of lower dimensional tori bend in a complex fashion, giving rise to the many tight loops seen in the picture. These manifolds are important because they are the objects able to carry the motion arriving along one of the resonances either to the "other part'' of the resonance or to a different resonance.

As long as the perturbation is increased, resonances are seen to become wider. Figure 2 displays the actual structure of action space (e-space) at two different (rather large) values of the perturbation parameter, namely  $\epsilon =0.02$ and  $\epsilon =0.04$. There, some details of the dynamics at moderate-to-high-level perturbations are shown in a plot similar to that in Fig. 1c; the character of the motion (resonant, quasi-periodic and stochastic) is represented in gray scale, from white to black, the different $\overline {Y}(t_{\rm f})$ intervals have been selected so as to sharpen the details of the phase space structure in each case.

Notice that while for  $\epsilon =0.02$ a significant part of the energy surface still looks regular, with wide resonance domains and a broad chaotic strip (Fig. 2a), it looks very chaotic for a somewhat slightly larger perturbation (Fig. 2b). Note however that within the weaker chaotic region, the MEGNO is still capable of unveiling the relics of resonance structures.

	Figure 3: Fraction of chaotic motion vs. $\epsilon$ (in logarithmic scale). The vertical lines refer to the values of $\epsilon$ adopted for Figs. 1c and 2.
Open with DEXTER

In Fig. 3 we present the fraction of chaotic motion as $\epsilon$ is increased, fixing a threshold of $\overline{Y}(t_{\rm f})=2$ below which orbits are regarded as regular. The global amount of stochasticity as a function of $\epsilon$is measured by counting how many pixels have a value of  $\overline {Y}$that exceeds the adopted threshold. The fraction of chaotic motion for the values of $\epsilon$ considered above, namely 0.02 and 0.04, are 38% and 91.7%, respectively. Note that there is some value of $\epsilon$, close to 0.03, where a transition to, let us say, global stochasticity takes place, in the sense that chaotic motion prevails in phase space, the chaotic zones not being necessarily connected. Further, for values of  $\epsilon\ge 0.1$ the system exhibits completely chaotic dynamics.

A quite remarkable result obtained by the application of the MEGNO to this model is its capability to display the actual dynamics at resonance crossings. We find that diffusion over the energy surface may occur, which is obvious for large perturbations, but in the case of smaller perturbations, the MEGNO shows a complex structure at the resonance intersections, which could restrain the spreading of chaotic motion, at least for the timescales considered.

3 Diffusion on the energy surface at moderate-to-high perturbations

We proceed to analyze the diffusion on the energy surface at moderate-to-high perturbations, following several orbits with initial conditions in different stochastic domains.

Table 1: Initial conditions on the energy surface for seven orbits with high values of the MEGNO. The energies h₁, h₂ are given in units of h/250, where $h=1/(4\beta ^4)$ is the total energy. The exact values of e₁ and e₂ should be obtained by means of Eqs. (8); their approximated values have been included so that they can be easily visualized on the energy surfaces.

For seven different initial conditions listed in Table 1 we have followed the wandering of the unperturbed integrals over the (e₁,e₂)-plane. We have picked two initial conditions in the chaotic strip close to h₁=0, (i) and (ii) - see Fig. 2a, another in the region surrounding the more regular central zone at the origin (iii) and four other initial conditions located near the crossing of resonances  (iv)-(vii); all of them have very high values of the MEGNO. The origin has also been considered; it actually corresponds to a regular orbit.

	Figure 4: Diffusion on the energy surface at moderate-to-high perturbations after $3\times 10^6$ characteristic periods T of the system for:   a)  $\epsilon =0.02$, and  b)  $\epsilon =0.04$; and after $3\times 10^8 T$ for:  c)  $\epsilon =0.02$, and  d)  $\epsilon =0.04$.
Open with DEXTER

The equations of motion have been integrated by means of a Runge-Kutta 7/8th order integrator (the so-called Dopri8 routine, see Hairer et al. 1987; Prince & Dormand 1981). The total motion time has been partitioned into 10⁵ intervals and the average of the unperturbed integrals over each interval has been obtained, with the aim of lowering the influence of bounded oscillations and so retaining the accumulated changes.

Let us remark that the integrator used for the present experiments is barely satisfactory since the resulting drift in the energy conservation is not negligible. However, the results obtained by means of such a standard routine are globally equivalent to those produced by the application of a rather more precise integrator, and succeed in providing a clear picture of the diffusion process in phase space over timescales like those considered here. We have already performed several runs using a 30th order Taylor scheme, which dramatically reduces the energy drift and is thus suitable for a long term diffusion investigation when integrations over much longer motion times are required (see Simó et al. 2004).

The results corresponding to $\epsilon =0.02$for a total motion time of  $3\times 10^6$ characteristic periods T of the system are presented in Fig. 4a (we refer to Fig. 2a for comparison). We note that even in the case of moderate perturbation, diffusion is completely irrelevant over such a timescale. Only when the resonances labeled by the harmonics (2,-1,0) and (2,0,-1) overlap does fast diffusion occur along such resonances, but it is restricted to a relatively small region of the energy surface. For the remaining initial conditions considered the unperturbed integrals remain confined to rather small domains. In particular, for the initial condition (vi) selected near the boundary h₃=0, the variation of the unperturbed integrals proves to be small, differing slightly from the expected behavior in case of stability. For the regular orbit at the origin, considered just for illustrative purposes, the wandering is restrained to a point, as expected.

The plot in Fig. 4b corresponds to $\epsilon =0.04$ and should be compared with the contour plot in Fig. 2b, where the chaotic regime prevails. Notice must be taken however that in the central region in the latter figure, for which  $t_{\rm f}=3500T$, there lie several chaotic orbits; the light gray points and even some of the white points within it have values of the MEGNO larger than 2 (as indicated in the caption). Recall that for constructing Fig. 2b the values of the MEGNO are binned in intervals chosen so as to highlight the dynamical structure of phase space, rather than discriminate stable quasiperiodic motion from chaotic motion. Thus, those chaotic orbits depicted in white in the contour plot having MEGNO values smaller than 2.1 are likely to behave in a regular fashion for rather short timescales, such as 3500T. Furthermore, Fig. 3 reveals that for  $\epsilon =0.04$, the system shows up as globally chaotic, with just a small fraction of phase space occupied by stable motion. This explains why, for the larger total motion time now considered, $3\times 10^6$ characteristic periods T of the system, diffusion manages to trespass the central zone, the points being mainly concentrated along a fairly well-defined strip corresponding to the curve associated to the (-2,1,1) resonance (see below). Note that the upper part of the central zone remains unvisited at this timescale, except for a rather strong but confined concentration of points near the corner.

On considering an even larger timescale,  $3\times 10^8$ characteristic periods T of the system, for $\epsilon =0.02$ (Fig. 4c) the fast diffusion occurring at the overlap of the resonances (2,-1,0) and (2,0,-1) after  $3\times 10^6$ periods now spreads and continues upwards along the resonances near the borders h₃=0 and h₂=0. For the initial condition (vi) chosen at a resonance crossing near h₂=0, the variation of the unperturbed integrals that was restrained to a small domain in Fig. 4a for the lower timescale, now moves upwards, then proceeds to the right and finally goes downwards, the path being drawn by three different resonances. Notice, however, that on this timescale the wandering does not reach the banana-shaped domain on the right corresponding to the initial condition (iii), probably due to the complexity of the dynamics at resonance intersections (recall how intricate such crossings may be). Summing up, on this timescale and for  $\epsilon =0.02$ the unperturbed integrals still roam over unconnected restraint zones of the energy surface. There remain very localized diffusion domains and the chaotic component seems far from being fully connected.

Instead, for $\epsilon =0.04$, after  $3\times 10^8$ characteristic periods T of the system the chaotic component is seen to be almost fully connected through the relics of the resonance structure (Fig. 4d). The upper right part of the energy surface, formerly empty, now appears densely populated by concentrations mainly along some still distinguishable remaining resonances.

Figure 5: a) Enlargement of Fig. 4d around (e1,e2)=(0,0). b) Diffusion of a single orbit over the energy surface for $\epsilon =0.04$ and $t_{\rm f}=6\times 10^8$ T (see text for details). c) Blow-up around the origin in the contour plot displaying the $\overline {Y}(t_{\rm f})$-levels on the energy surface for $\epsilon =5\times 10^{-3}$ and  $t_{\rm f}=350T$ (taken from Cincotta et al. 2003). d) Diffusion close to (e1,e2)=(0,0) for $\epsilon =5\times 10^{-3}$ and $t_{\rm f}=3\times 10^8$ T (see text for details).

Open with DEXTER

An enlargement of the central region in Fig. 4d is displayed in Fig. 5a, to see whether any unvisited domain still remains in that part of the energy surface for such a large timescale. A strong concentration of points persists along the strip already present in Fig. 4b for $3\times 10^6$ periods, namely, the curve associated to the (-2,1,1) resonance. Some further noticeable concentrations along an arc near the upper corner of the figure and several straight lines give a clear account of the remaining resonance structure.

In view of the results obtained for $\epsilon =0.04$ on the larger timescale, we integrate a single orbit, the initial condition (i) taken in the chaotic strip near the border h₁=0, for an even larger $t_{\rm f}$, $\;6\times 10^{8}$ characteristic periods T, and trace the path of its unperturbed integrals over the energy surface, with the aim of verifying whether this path covers the whole chaotic component. The results are displayed in Fig. 5b, and should be compared with the plot in Fig. 4d for seven different initial conditions in Table 1 corresponding to chaotic orbits. As seen in the picture, for this single orbit the unperturbed integrals trail over the chaotic component, the relics of resonances serving as routes. It should be mentioned that, for half the total motion time considered, the upper part of the energy surface remained yet unreached, and it took a much longer time to visit the whole chaotic component.

Finally, and in order to confirm the results presented in Fig. 5c - obtained for $\epsilon =5\times 10^{-3}$ and a rather short time motion $t_{\rm f}=350T$, and already discussed in Sect. 2 - in the sense that the intricate character of the dynamics at resonance crossings may well restrain diffusion, we have selected an initial condition near the origin (e₁,e₂)=(0.076342430317, -0.0799030662741) with a rather high value of the MEGNO, namely 5.89, and follow its path for over $3\times 10^8$ T (see Fig. 5d). We observe that despite the very large timescale considered, the variation of the unperturbed integrals remains confined to a rather small domain on the energy surface.

Note that when dealing with chaotic systems the use of different processors and/or software may lead to somewhat different results, though the global structures obtained on the energy surface would not show significant differences. It might happen, however, that two very close chaotic domains could show up on the energy surface either connected or unconnected when the integrations are carried out with different processors.

  
4 Discussion

In the present paper we present a preliminary study of diffusion in phase space for a rather simple dynamical system in order to elucidate its efficiency in connecting different chaotic components. A thorough investigation of long term diffusion as well as of the efficiency of a single orbit in exploring the whole energy surface in the case of a fully chaotic dynamics was provided in Simó et al. (2004), where particular features of the dynamics of the chosen system are included.

Despite the simplicity of the adopted model, several results concerning its dynamics apply to any 3D Hamiltonian system exhibiting a divided phase space, such as galactic or planetary systems. In fact, two values for the perturbation have been selected such that the dynamics of our toy model resembles that of a galactic system (a moderate perturbation, in which case the two components are comparable) and asteroidal dynamics (a large perturbation for which the chaotic component prevails).

We have observed that diffusion may actually take place over the full chaotic component only for strongly chaotic dynamics and very long timescales. Notice that while timescales of the order of $\sim$ $6~\times~ 10^8$ characteristic periods are physically significant in asteroidal dynamics, they are completely irrelevant for galactic systems, for which $\sim$10³ periods is a realistic upper bound. Anyway, despite the motion time elapsed, for moderate perturbations diffusion seems to be unable to connect different stochastic domains.

Since galaxies are N-body systems, one should not expect to gain relevant information on the time evolution of the distribution by following a single star orbit over a long stellar time. But from the above results, a single connected chaotic domain should exist for galaxies having a very large chaotic component (similar to that given in Fig. 2b); then it could happen that the timescale for diffusion given above ($\sim$10⁸ periods) is highly overestimated, the actual timescale required for the galaxy to reach a nearly stationary distribution being somewhat lower.

An interesting fact is that diffusion is driven along the resonances, even in an overlap regime. Though the resonance structure is destroyed by the overlap, diffusion proceeds along its relics; the latter can be clearly seen in the MEGNO contour plots. This kind of diffusion should not be attributed strictly to Arnol'd diffusion (Arnol'd 1964; Chirikov 1979; Cincotta 2002), where the existence of a chain of tori (the so-called transition chain) is required to have a path for a single orbit to transit along the resonance. However, according to the obtained numerical evidence, the mechanism leading to diffusion shows a somewhat geometrical resemblance to Arnol'd diffusion in the sense that the strong stochastic layers around resonances serve as routes for diffusion, thus managing to connect two distant chaotic domains of phase space.

We have shown the complex nature of dynamics at resonance intersections in a multidimensional system. Moreover, we have provided numerical evidence that the predicted stability domains in the resonance crossings would restrain diffusion, at least over rather large motion times. Nonetheless, the evidence obtained here does not invalidate the theoretical conjecture that diffusion may spread over the whole energy surface of multidimensional systems through resonance intersections, which still remains an open question.

Observational evidence suggests that a model resembling for instance an elliptical galaxy should exhibit a phase space structure analogous to that presented in Fig. 2a (corresponding to a moderate perturbation). In such a case, we have observed that the unperturbed integrals remain confined to rather small chaotic domains and that diffusion is completely incapable of connecting them, at least over timescales of the order of $\sim$10⁸ characteristic periods, several orders of magnitude larger than the typical timescales for galactic systems. Further, the existence of remarkable point concentrations on the energy surface suggests that the distribution function is unlikely to depend only on the energy level, but also on one or two more local pseudo-invariants in a somewhat complicated fashion. This is a natural consequence of the discontinuous dependence of the unperturbed actions on the phase space coordinates in nearly integrable Hamiltonian systems.

This investigation encourages the study of more realistic models. In this direction, some very preliminary results concerning a perfect ellipsoid perturbed by Denhen's law for modeling an elliptical galaxy agree with those obtained here for the toy model. These issues will be thoroughly addressed in a forthcoming paper.

Acknowledgements

The authors would like to acknowledge the CONICET (Argentina) for financial support through the Instituto de Astrofísica de La Plata. We are also grateful to the valuable comments of C. Simó. Finally, we appreciate the referee labour on reviewing our manuscript in a rather critical fashion, thus contributing to its improvement.

Appendix A: 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 Cincotta et al. (2003). 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 ${\bm{p}},~{\bm{q}}\in {\rm I\!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 ${\bm{x}}=({\bm{p}},{\bm{q}})\in {\rm I\!R}^{2N}\hspace{-1mm},~ {\bm{v}}=(-\partial H/\partial{\bm{q}},\ \partial H/\partial {\bm{p}})\in{\rm I\!R}^{2N}$, the equations of motion read

$$\dot{\bm{x}}={\bm{v}}({\bm{x}}).$$ (A.1)

Let $\gamma(t)$be an arc of an orbit of flow (A.1) on a compact energy surface $M_h\subset{\rm I\!R}^{2N}$ , $M_h=\{{\bm{x}}: H({\bm{p}},{\bm{q}})=h\}$, so that

$$\gamma(t)=\left\{{\bm{x}}(t';{\bm{x}}_0):{\bm{x}}_0\in M_h,\ 0\le t'< t\right\},$$ (A.2)

and the full positive orbit is $\gamma=\lim_{t\to\infty}\gamma(t)$.

Relevant information about the flow in the vicinity of any orbit $\gamma$ is gained through its largest Lyapunov Characteristic Number (LCN-hereafter), $\sigma(\gamma)$, defined as

$$\sigma(\gamma)=\lim_{t\to\infty}\sigma_1(\gamma(t)),\quad \si... ...ert\bm{\delta}(\gamma(t))\Vert\over\ \Vert\bm{\delta}_0\Vert},$$ (A.3)

with $\bm{\delta}(\gamma(t))$ and $\bm{\delta}_0$ "infinitesimal displacements" from $\gamma$ at times t and 0, respectively, (see below) and where $\Vert\cdot\Vert$ is some norm. The fact that the LCN measures the "mean exponential rate of divergence of nearby orbits", is stated explicitly when recasting (A.3) in the integral form

$$\sigma(\gamma)=\lim_{t\to\infty}{1\over t}\int_0^t{{\dot{\del... ...(t'))}{\rm d}t'}= \overline{\left(\dot{\delta}/\delta\right)},$$ (A.4)

with $\delta\equiv\Vert\bm{\delta}\Vert,~\dot{\delta}\equiv {\rm d}\delta/{\rm d}t= \dot{\bm{\delta}}\cdot\bm{\delta}/\Vert\bm{\delta}\Vert$, the bar denoting time-average. Recall that the tangent vector  $\bm{\delta}$ satisfies the variational equation

$$\dot{\bm{\delta}}=\Lambda(\gamma(t))\bm{\delta},$$ (A.5)

where $\Lambda=D\bm{v}$ is the Jacobian matrix of the vector field v.

We now introduce the MEGNO, $Y(\gamma(t))$, through the expression

$$Y(\gamma(t))={2\over t}\int_0^t{\dot{\delta}(\gamma(t'))\over\delta(\gamma(t'))}t'{\rm d}t',$$ (A.6)

which is related to the integral appearing in (A.4). Notice that in the case of an exponential increase of $\delta$, $\delta(\gamma(t))=\delta_0 \exp(\lambda t)$, the quantity  $Y(\gamma(t))$ can be considered as a weighted variant of the integral in (A.4). Instead of the instantaneous rate of increase, $\lambda$, we average the logarithm of the growth factor, $\ln~(\delta(\gamma(t))/\delta_0)=\lambda 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 Cincotta et al. (2003), for any such (quasi-periodic) orbit, $\gamma_{\rm q}$, the temporal evolution of $Y(\gamma_{\rm q}(t))$ is given by

$$Y\left(\gamma_{\rm q}(t)\right)\approx 2-{\ln~(1+\lambda_{\rm... ...^{2}\over\lambda_{\rm q} ~t}+ O\left(\gamma_{\rm q}(t)\right),$$ (A.7)

where O denotes an oscillating term with zero average. The $\lim_{t\to\infty}Y\left(\gamma_{\rm q}(t)\right)$ does not exist but, on introducing the timeaverage

$$\overline{Y}(\gamma_{\rm q}(t))\equiv{1\over t}\int_0^tY(\gamma_{\rm q}(t')){\rm d}t',$$ (A.8)

it can readily be shown that

$$\overline{Y}\left(\gamma_{\rm q}\right)\equiv\lim_{t\to\infty}\overline{Y} \left(\gamma_{\rm q}(t)\right)=2.$$ (A.9)

Therefore, for the case of quasi-periodic motion, $\overline{Y}(\gamma)$ is a fixed constant, independent of $\gamma$.

For an irregular orbit, $\gamma_{\rm i}$, within any stochastic component, for which $\delta(\gamma_{\rm i}(t))\approx\delta_0 \exp{(\sigma_{\rm i}t)}$ where $\sigma_{\rm i}$ is $\gamma_{\rm i}$'s LCN, the temporal evolution of the MEGNO reduces to

$$Y\left(\gamma_{\rm i}(t)\right)\approx\sigma_{\rm i}t+\tilde{O}\left(\gamma_{\rm i}(t)\right),$$ (A.10)

where $\tilde{O}$ 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

$$\overline{Y}\left(\gamma_{\rm i}(t)\right)\approx{\sigma_{\rm i}\over 2}~t,\quad t\to \infty.$$ (A.11)

Therefore, for a chaotic orbit both $Y\left(\gamma_{\rm i}(t)\right)$ and $\overline{Y}\left(\gamma_{\rm i}(t)\right)$ 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 $\delta$ grows with some power of t, say n, where $\overline{Y}\to 2n$ as  $t \to \infty$.

Therefore, the MEGNO's temporal evolution allows of being summed up in a single expression valid for any kind ofmotion. The asymptotic behavior of  $\overline{Y}(\gamma(t))$may be written in the fashion  $\overline{Y}(\gamma(t))\approx a_{\gamma}t+d_{\gamma}$, where  $a_{\gamma}=\sigma_{\gamma}/2$ and  $d_{\gamma}\approx 0$ for irregular, stochastic motion, while  $a_{\gamma}=0$ and  $d_{\gamma}\approx 2$ for stable quasi-periodic motion. Departures from the value  $d_{\gamma}\approx 2$ indicate that $\gamma$ is close to some periodic orbit, being $d_{\gamma}\la 2$and  $d_{\gamma}\ga 2$ for stable and unstable periodic orbits, respectively (see Cincotta et al. 2003 for details). Note also that the quantity $\hat{\sigma}_1=Y/t$ verifies

$$\hat{\sigma}_1(\gamma_{\rm q}(t))\approx{2\over t},\qquad\hat... ..._1(\gamma_{\rm i}(t)) \approx\sigma_{\rm i},\qquad t\to\infty,$$ (A.12)

which shows that in the case of regular motion $\hat{\sigma}_1$ converges to 0 faster than $\sigma_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  $\overline {Y}$ are given in Cincotta & Simó (2000) and Cincotta & Giordano (2002). An interesting application to a 2.5D problem, namely, the Arnold's classical problem of diffusion when the two small parameters are equal, may be found in Simó (2001). The associated splitting of separatrices has been studied in Simó & Valls (2001). Recent applications of the MEGNO to Solar System dynamics may be found in, for instance Gozdziewski (2002, 2003a,b). For a generalized version of the MEGNO and its application to a 2D and a 4D area preserving maps we refer to Cincotta et al. (2003).

References

• Arnol'd, V. I. 1964, Sov. Math.-Dokl., 5, 581 [NASA ADS] (In the text)
• Cincotta, P. M. 2002, New Astron. Rev., 46, 13 [NASA ADS] (In the text)
• Cincotta, P. M., & Giordano, C. M. 2002, in Advances in space dynamics 2. Applications in Astronomy, ed. O. Winter, & A. F. B. A. Prado, 237 (In the text)
• Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Phys. D, 182, 151 [NASA ADS] (In the text)
• Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [NASA ADS] (In the text)
• Chirikov, B. V. 1979, PhR, 52, 263 [NASA ADS] (In the text)
• Gozdziewski, K. 2002, A&A, 393, 997 [EDP Sciences] [NASA ADS] (In the text)
• Gozdziewski, K. 2003a, A&A, 398, 315 [EDP Sciences] [NASA ADS] (In the text)
• Gozdziewski, K. 2003b, A&A, 398, 1151 [EDP Sciences] [NASA ADS]
• Hairer, E., Nørsett, S., & Wanner, G. 1987, Solving ordinary differential equations I: Nonstiff problems (New York: Springer-Verlag) (In the text)
• Merritt, D. 1999, PASP, 111, 129 [NASA ADS] (In the text)
• Merritt, D., & Fridman, T. 1996, ApJ, 460, 136 [NASA ADS] (In the text)
• Merritt, D., & Valluri, M. 1996, ApJ, 471, 82 [NASA ADS] (In the text)
• Prince, P., & Dormand, J. 1981, J. Comp. Appl. Math., 35, 67 [NASA ADS] (In the text)
• Simó, C. 2001, in Global analysis of dynamical systems, ed. H. W. Broer, B. Krauskopf, & G. Vegter (IOP Publishing), 263 (In the text)
• Simó, C., Giordano, C. M., & Cincotta, P. M. 2004, in preparation (In the text)
• Simó, C., & Valls, C. 2001, Nonlinearity, 14, 1707 [NASA ADS] (In the text)
• Valluri, M., & Merritt, D. 1998, ApJ, 506, 686 [NASA ADS]