4 Measure of chaos and fine structure of the phase space

  
4.1 The Mean Exponential Growth factor of Nearby Orbits (MEGNO)

In Sect. 2 we have shown that if the potential is not too far from the central symmetry, we can derive the global orbital structure by means of an additional pseudo-invariant. With this model we cannot describe chaotic zones (like the stochastic layer that separates box and loop families) and secondary resonances (like the banana sub-family). Hence to get information about the fine structure of the phase space as well as to measure chaos in the irregular components we should follow some other approach. As we pointed out in Sect. 1 the LCN provides a good measure of chaos but it does not furnish any information about the structure of the regular component. This fact is a consequence of the definition of the LCN, $\sigma$, for a given orbit, $\gamma(t)$, on a compact energy surface, $M_{\rm h}$,

$$% \sigma(\gamma)=\lim_{t\to\infty}{1\over t}\ln{\delta\left(\gamma(t)\right) \over\delta_0},$$ (21)

(if the limit exists) where $\delta(\gamma(t))\equiv\vert\vec{\delta}(\gamma(t))\vert$ and $\delta_0$ are infinitesimal displacements from $\gamma$ at times t and 0, respectively. $\vec{\delta}(\gamma(t))$ satisfies the variational equations $\dot{\vec{\delta}}=\Lambda(\gamma(t))\,\vec{\delta}$, where $\Lambda$ is the Jacobian matrix of the vector field.

In practice, instead of the infinite limit one computes $\sigma(\gamma(T))$ for T "large enough'', for instance, 10⁴ characteristic periods. For any regular orbit, $\delta(\gamma(t))\approx\delta_0(1+\lambda_{\gamma}t)$ - see below - and then $\sigma(\gamma(T))\approx\ln T/T$, almost independent of $\gamma$. However some departures from the latter value exist depending whether the orbit is periodic stable or it is stable quasiperiodic but passing close to an unstable periodic orbit.

As an example let us consider a neighbourhood, U, of an unstable periodic orbit $\gamma_{\rm u}$. The motion in U is mainly determined by the stable and unstable manifolds associated to $\gamma_{\rm u}$. Therefore any quasiperiodic orbit, $\gamma_{\rm q}$, that falls in U will mimic, for a short time interval, $\gamma_{\rm u}$. Since for the latter $\delta(\gamma_{\rm u}(t))$ grows exponentially with time, then $\delta(\gamma_{\rm q}(t))$ will behave in a similar manner while $\gamma_{\rm q}$ lies in U. This will happen periodically, each time $\gamma_{\rm q}$enters in U. Therefore if we look at the time evolution of $\sigma(\gamma_{\rm q}(t))$ we should see a $\ln t/t$ law modulated by a periodic peak structure (besides the purely quasiperiodic oscillations). Anyway the final value of $\sigma(\gamma_{\rm q}(t))$ for t=T will be very close to $\ln T/T$, except (perhaps) if we stop the computation just when $\gamma_{\rm q}$ is in U. The only way to put in evidence the latter behaviour is to take into account somehow the time evolution of $\sigma(\gamma(t))$. A simple way to do that and to amplify the effect of unstable periodic orbits on quasiperiodic motion is the following.

The definition of the LCN, given by (21), can be written in an integral form as

$$\begin{eqnarray*}\sigma(\gamma)=\lim_{T\to\infty}{1\over T}\int_0^T{\dot{\delta}... ...\right)}\,{\rm d}t\equiv \left<{\dot{\delta}\over\delta}\right>, \end{eqnarray*}$$

where $\dot{\delta}\equiv\dot{\vec{\delta}}\cdot\vec{\delta}/\delta$ is the time derivative of $\delta(\gamma(t))$ and $\langle \cdot \rangle$ denotes the usual time-average. Let us define the Mean Exponential Growth factor of Nearby Orbits, ${\cal J}$, as

$$% {\cal J}(\gamma(T))\equiv{2\over T}\int_0^T{\dot{\delta}\le... ...amma(t)\right) \over \delta\left(\gamma(t)\right)}\,t{\rm d}t.$$ (22)

Then for any regular, quasiperiodic orbit, $\gamma_{\rm r}$, we get

$$% {\cal J}(\gamma_{\rm r}(T))\approx 2\left(1-{\ln(1+\lambda_... ...ambda_{\gamma_{\rm r}}T}\right)+ {\rm O}\,(\gamma_{\rm r}(T)),$$ (23)

where, as before, $\lambda_{\gamma_{\rm r}}>0$ is the linear rate of divergence in a neighbourhood of $\gamma_{\rm r}$ and ${\rm O}\,(\gamma_{\rm r}(T))$ denotes an oscillating term $(\langle {\rm O}\,(\gamma_{\rm r}(T))\rangle=0)$ introduced by the corresponding quasiperiodic terms in $\delta(\gamma_{\rm r}(T))$. Indeed, for stable, regular motion $\delta(\gamma(t))$ can be approximated by $\delta\approx\delta_0[1+\lambda t+ t\,u(t)]$ where u is an oscillating function of t and |u(t)| is bounded, say $\vert u(t)\vert< b<\lambda$. So it is straightforward to show that for large T, ${\cal J} (\gamma_{\rm r}(T))$ oscillates about 2 and when $T\to\infty$ with an amplitude

$$\begin{eqnarray*}\vert{\cal J}(\gamma_{\rm r}(T))-2\vert<4\ln{\lambda_{\gamma_{\rm r}}+b\over \lambda_{\gamma_{\rm r}}-b}\cdot \end{eqnarray*}$$

Hence if the motion is quasiperiodic, then ${\rm O}(\gamma_{\rm r}(T))$ represents eventually small oscillations around the mean value of ${\cal J}(\gamma_{\rm r})\equiv {\cal J}_{\rm r}$, that can be neglected. In case $\gamma_{\rm r}$ is close to a periodic orbit, ${\rm O}(\gamma_{\rm r}(T))$ exhibits a nearly periodic character and small amplitude.

As we note from (23), $\lambda_{\gamma_{\rm r}}$ determines how fast ${\cal J}_{\rm r}\rightarrow 2$. Since $\lambda$ is a measure of the lack of isochronicity around $\gamma$ (to be precise, $\lambda$ is the largest eigenvalue of the matrix $\partial\vec{\omega}/\partial\vec{I}$, being $\vec{I}$the action), we see that the smaller $\lambda$ appears for $\gamma$ in a neighbourhood of a stable periodic orbit. Therefore we expect a slower rate of convergence of ${\cal J}_{\rm r}$ to 2 for $\gamma_{\rm r}$ close to a stable "elliptic'' periodic orbit. On the other hand, from the discussion given above about the behaviour of quasiperiodic orbits close to unstable periodic ones, we expect that in this case ${\cal J}_{\rm r}$ presents quasiperiodic oscillations, given by ${\rm O}(\gamma_{\rm r}(T))$, as well as a peak structure. Therefore for arbitrary $\gamma_{\rm r}$, the formal limit, $\lim_{T\to\infty}{\cal J}(\gamma_{\rm r}(T))$, does not exist (see below).

For any irregular orbit, $\gamma_{\rm i}$, for which $\delta$ grows exponentially with time, we get

$$% {\cal J}_{\rm i}\equiv{\cal J}(\gamma_{\rm i}(T))\approx\sigma_{\gamma_{\rm i}} T+ \widehat{{\rm O}}(\gamma_{\rm i}(T)),$$ (24)

where $\sigma_{\gamma_{\rm i}}$ is the LCN for $\gamma_{\rm i}$ and $\widehat{{\rm O}}(\gamma_{\rm i}(T))$ is an oscillating (in general neither quasiperiodic nor periodic) term.

While ${\cal J}$ has not a formal limit, ${\cal J}/T$ and the mean values $\bar{{\cal J}}_{\rm r}, \bar{{\cal J}}_{\rm i}$ have an asymptotic law for $T\rightarrow\infty$

$$% \bar{{\cal J}}_{\rm r}\approx 2,\quad \bar{{\cal J}}_{\rm i... ...\over T},\quad {{\cal J}_{\rm i}\over T}\approx\sigma_{\rm i}.$$ (25)

For regular motion $\bar{{\cal J}}_{\rm r}$ is constant, almost independent of initial conditions and any other parameter while, for irregular motion, $\bar{{\cal J}}_{\rm i}$ grows linearly with time with a rate that is the LCN. Note that $\hat{\sigma}_{\rm r}={\cal J}_{\rm r}/T$ converges faster to zero than $\sigma_{\rm r}(T)= \ln T/T$, while $\hat{\sigma}_{\rm i}={\cal J}_{\rm i}/T$ converges to the actual non-zero ${\rm LCN}$, $\sigma_{\rm i}$, with the same rate as it does $\sigma_{\rm i}(T)$(compare Eqs. (24), (25) with (21)).

Of most importance is that $\bar{{\cal J}}(T)$ can be written in an unique way for both types of motion, $\bar{{\cal J}}(T)\approx aT+b$ with a=0, b=2 for regular, quasiperiodic motion and $a=\sigma$, $b=b_0\approx 0$ for irregular one. If $\delta(\gamma(t))$ grows with some power of t, say n, as it could happen in some degenerated cases, a=0, b=2n. Only when the phase space has a hyperbolic structure, where nearby orbits diverge exponentially with time, $a\ne 0$ and the MEGNO grows with time. This occurs for irregular, chaotic motion and also, for instance, for unstable periodic orbits (see Giorgilli et al. 1997, for a visualization of the hyperbolic structure of chaotic zones).

From (25) it turns out that if we have ${\cal J}(T)$ for any T we can recover the LCN by a linear least squares fit. The main advantage of this approach is that we use the dynamical information contained in ${\cal J}(T)$ along the whole time interval. Hence, we expect that this procedure will provide a good estimation of the LCN in both regular and irregular domains. Furthermore a least squares fit will also give us information about the location of hyperbolic orbits, the very origin of chaos. Therefore the derivation of the LCN from the MEGNO seems to be useful not only to get the global dynamics but to learn some details concerning the fine structure of phase space as well. A more complete discussion about the MEGNO technique as well as some examples of application to 3D systems are given in Cincotta et al. (2000).

  
4.2 Application to the 2D logarithmic potential

Figure 7: Time evolution of $\bar{{\cal J}}$ for regular and irregular orbits. On the left, for two regular, stable quasiperiodic orbits, A, B, but orbit B is very close to a unstable periodic orbit. On the right, we plot both $\bar{{\cal J}}$ (smooth curve - actually $2\bar{\cal J}$, see text) and ${\cal J}$ (noisy curve) for an irregular orbit, C, in the stochastic layer. All these orbits correspond to the logarithmic Hamiltonian (16) with p02=2, $r_{\rm c}=0.1$, q=0.9, and on the energy surface $M_{\rm h}$ defined by h=-0.4059 and initial conditions $x_0\approx 0.33$, 0.02, 0.002, y0=0, px=0 for A, B and C, respectively

To investigate numerically this technique we considered again the logarithmic potential defined in (15), or its associated Hamiltonian (16), for the same parameters, energy level, $M_{\rm h}$, and q values taken in Sect. 2.2. We use the initial values x₀, p_x0 to parameterize loop and box orbits respectively (see Sect. 2.2 and below).

The computation of ${\cal J}$ was done using (22) for a given set of initial conditions. All the integrations were carried out for a realistic time scale, $T\sim 10^3\ T_{\rm D}\approx 3000$ where $T_{\rm D}$ is the period of the long-axis periodic orbit. Therefore the computational effort per unit time is almost the same needed to compute the LCN but comparatively shorter motion times are required. The renormalization of $\vec{\delta}$ (if necessary), proceeds naturally from (22).

To solve the variational equations we took $\vec{\delta}_0$ along the x axis for loops and along the p_x axis for boxes with $\vert\vec{\delta}_0\vert=1$ and random sign. We used a Runge-Kutta 7/8 th order integrator (the so-called DOPRI8 routine, see Prince & Dormand 1981; Hairer et al. 1987). The accuracy in the preservation of the value of the energy is $\sim 10^{-13}$.

To eliminate fast quasiperiodic oscillations, that is, to compute $\bar{{\cal J}}(T)$, we averaged ${\cal J}(T)$ as follows

$$% \bar{\cal J}(T_k)={1\over k}\sum_{i=1}^k{\cal J}(T_{i}),\quad T_k=T_0+k\Delta T,$$ (26)

where $\Delta T\approx 0.06$ is the time-step. Hence $\bar{{\cal J}}$ depends on T and weakly on $\Delta T$. Alternatively, if necessary, $\bar{{\cal J}}(T)$ computed as $T^{-1}\int_0^T {\cal J}(t){\rm d}t$ would provide a smoother behaviour than (26), which is independent of the time step.

To perform the least squares fit to get the slope of ${\cal J}(T)$, that is the LCN, we use the values of $\bar{{\cal J}}(T)$ along the last 85% of the time interval ( $450\le T\le 3000$), just to avoid the initial transient. We add a factor 2 in the derived slope to compensate the average introduced in ${\cal J}$. Indeed, since for an irregular orbit ${\cal J}$ grows nearly linear, the slope derived from $\bar{{\cal J}}$ would be underestimated in a factor 2.

In Fig. 7 we show the time evolution of $\bar{{\cal J}}$ for three representative orbits. The regular ones belong to the loop family while the irregular one to the stochastic layer. For the latter we plot ${\cal J}$ together with $2\bar{\cal J}$ to put in evidence that the factor 2 introduced ad hoc is necessary. The figure on the left corresponds to two regular orbits, A and B. Orbit A, stable quasiperiodic, saturates very fast from below to 2, without any significant oscillation. Orbit B, also stable quasiperiodic, comes very close to an unstable periodic orbit. We observe the influence of the unstable periodic orbit on B leading to several local maxima of decreasing amplitude. This behaviour of the amplitude of the maxima is due to the average of ${\cal J}$ (see Eq. (26)) and it should decrease as $\sim 1/T$. Note that in this case $\bar{{\cal J}}$ takes higher values. This is also due to the presence of a nearby hyperbolic orbit. For the irregular orbit, C, (on the right figure) we see a nearly linear behaviour. In fact, $\bar{{\cal J}}(T)$ looks smoother and follows the same linear trend than ${\cal J}(T)$.

In Fig. 8 we plot the time evolution of $\hat{\sigma}(T) \equiv{\cal J}/T$ and $\sigma(T)$ for orbits A and C. For the regular orbit A, $\hat{\sigma}(T)<\sigma(T)$. The theoretical values (in logarithmic scale) are -3.18 and -2.57 respectively (see Eq. (25) and around), which are in good agreement with the computed ones. For the irregular orbit C, we see that $\hat{\sigma}(T)\approx\sigma(T)$, also consistent with the above discussion.

Figure 8: Time evolution of $\log\hat{\sigma}(T)$, $\log\sigma(T)$ for orbit A, stable quasiperiodic and C, irregular, shown in Fig. 7 (see text)

Since the relative error in the estimation of the positive LCN after a motion time T is $\sim T_{\rm L}/T$, where $T_{\rm L}=1/\sigma$ is the Lyapunov time, we see that $T\sim 10^3\ T_{\rm D}$ is not enough to separate a chaotic region with $T_{\rm L}\sim 10^3\ T_{\rm D}$ from the regular one. This is one of the reasons of why it is necessary to take very long motion times to compute the LCN using the standard method. It is important to remark that in this particular application it is enough to obtain an estimation of the order of magnitude of $T_{\rm L}$. When an accurate determination of the LCN is necessary, the motion time could be very large. The presence of small resonance domains embedded in a chaotic sea produces the so-called stickiness that reduces the free diffusion. So the motion time needed in this case to compute the positive LCN should be large enough so that the orbit could fill almost all the available subset of the energy surface.

Figure 9: $\log\sigma$, $\log\sigma_{\rm ls}$ for the domains of loop orbits (left column) and box orbits (right column). Each point represent an orbit labeled by the initial value of x-coordinate, x0, for loops and the initial value of px-coordinate, px0, for boxes. The ensembles include, in each case, about 3500 orbits and the total motion time is $T=3000\approx 10^3\ T_{\rm D}$ (see text for details)

Let us now consider ensembles of orbits. To explore the phase space we use a 1D initial conditions space. The best choice would be, for instance, points along a maximum circle connecting the NP and SP of the sphere shown in Fig. 6 and passing through the 1:1 and 2:1 periodic orbits. In this way we would include loops, boxes, bananas, etc. Thus, the angle $\varpi$ is a good parameter to label orbits in the main families. Nevertheless, to compare our results with that given by PL96 we follow their approach. Thus, for loops we take the ensemble $L_{\rm q}=\{0\le x_0\le X(q),\, y_0=0,\, p_{x_0}=0\}\subset M_{\rm h}$where X(q) corresponds to the location along the x-axis of the 1:1 periodic orbit for a given value of q. In the same way, for box orbits, $B_{\rm q}=\{x_0=0,\, y_0=0,\, 0\le p_{x_0}\le P(q)\}\subset M_{\rm h}$where $P(q)\equiv P$ corresponds to the location along the p_x-axis of the long-axis periodic orbit. For the values of q considered here, one numerically finds that $X(q)\approx 0.44$, 0.38, 0.33 for q=0.9, 0.8, 0.7, respectively (compare with the values derived in Sect. 2.2) while $P=(2h-p_0^2\ln r_{\rm c}^2)^{1/2}\approx 2.9$. We take ensembles of about 3500 initial conditions, similar to that considered by PL96. For each orbit we compute ${\cal J}(T)$, $\bar{{\cal J}}(T)$, $\sigma(T)$ after $T=3000\approx 10^3\ T_{\rm D}$ and we derive the slope of $\bar{{\cal J}}(T)$, $\sigma_{\rm ls}$, by a least squares fit in the way described above.

In Fig. 9 we show in the same plot $\sigma$ and $\sigma_{\rm ls}$ for $L_{\rm q}$ and $B_{\rm q}$ (q=0.9, 0.8 and 0.7). First of all we note that in any case $\sigma$ and $\sigma_{\rm ls}$ agree in the gross stochastic layer. The same happens, for q=0.7, in the thin stochastic layer around the x-axis periodic orbit as well as in some narrow chaotic zones around other resonances. But the value of $\sigma$ over all regular regions is nearly the same, $\log\sigma\approx\log(\ln T/T)\approx-2.57$. Just a few zones, where $\sigma$ appears to behave in a smoother way, could be supposed to be the signature of an island. In contrast, $\sigma_{\rm ls}$ clearly shows the underlying structure of the regular region. Note that $\sigma_{\rm ls}$ leads to a Lyapunov time for the regular component of $T_{\rm L}\sim 10^6\ T_{\rm D}$ in $10^3\ T_{\rm D}$ while $\sigma$ leads to $T_{\rm L}\sim 10^2\ T_{\rm D}$. To get such long values of $T_{\rm L}$ (in the regular domain) by means of the computation of $\sigma$, the total motion time should be $T\sim 10^7\ T_{\rm D}$.

Figure 10: $\log\sigma_{\rm ls}$ for the region of loop domain occupied by resonances. All significant resonances are marked by a full-line arrow. Those indicated by a dashed-line arrow and labeled by "r'' are also detected by PL96 by means of the FMA (see text)

Globally we see that, for $L_{\rm q}$ (Fig. 9 - left column), the domain is clearly divided in two zones. One near to the unstable short-axis orbit (at the origin), that contains irregular orbits and many small sub-families corresponding to each small resonance domain (see below). The other zone, near to the 1:1 periodic orbit, looks free of resonances and completely populated by quasiperiodic loop orbits. Note that $\sigma_{\rm ls}$ increases slowly as we approach to the 1:1 orbit. This is a consequence of the fact that the rate of convergence of ${\cal J}$ to 2 is the slower for orbits near to the latter "elliptic'' periodic orbit (see Sect. 4.1).

Figure 11: a) Section y=0 for loops close to the border of the stochastic layer, window: $0.04\le x \le 0.076$, $\vert p_x\vert\le 0.08$. The arrows indicate some of the resonances marked in Fig. 10 (q=0.7) near to the stochastic layer. b) $\bar{{\cal J}}$ after T=3000 and for x0 in the same interval as above. The range in $\bar{{\cal J}}$ is [1.98,2.035] and the dotted line is the level $\bar{\cal J}=2$ (see text)

For $B_{\rm q}$ (right column) the scenario is different. The fraction of irregular orbits is larger (specially for q not too far from 1) and we can appreciate many resonances along the whole domain, even for q=0.9. For the largest perturbation, q=0.7, one could infer that quasiperiodic box orbits do exist but, in some sense, discontinuously and the region occupied by other sub-families (boxlets) is almost as large as that filled by purely box orbits. Recall that the banana sub-family is not included here because all that orbits cross the surface y=0 near the boundary of the section with $x\ne 0$ (this follows immediately from Fig. 1; see also Fig. 15). As a difference with $L_{\rm q}$, $\sigma_{\rm ls}$ decreases slowly as we approach to the long-axis periodic orbit. This behaviour is purely due to the choice of variables, since if instead of (x,p_x) we had used (y,p_y), this orbit would appeared as an elliptic point. This effect is a consequence of the projection of the motion onto a 2D plane while, if we project it onto the "natural'' space, which is the 2D sphere, this problem does not appear.

Figure 12: $\log\sigma_{\rm ls}$ for the whole box domain but in different windows. All significant resonances are marked by a full-line arrow. Those indicated by a dashed-line arrow and labeled by "r'' are also detected by PL96 by means of the FMA (see text)

In Fig. 10 we plot the resonance zone for loops close to the stochastic layer. We marked by an arrow the most significant resonances. Those labeled by "r'' were also detected by PL96 by means of the frequency map analysis (FMA) and therefore we have at hand the rotation number corresponding to these resonances (see Sect. 6). It is important to mention that the resolution in x₀ is similar to that of PL96 but the motion time used here is larger, about a factor 10 ( $10^3\ T_{\rm D}$ and $10^2\ T_{\rm D}$, respectively). Yet, Fig. 7 shows that shorter motion times, about $500\ T_{\rm D}$, would be enough. Note that, in any case, this technique is able to put in evidence the complex structure of resonances in the neighbourhood of the stochastic layer. We distinguish basically two different types of departures from quasiperiodic motion, peaks and valleys. From the discussion given in Sect. 4.1 we easily conclude that the peaks correspond to unstable periodic orbits and to quasiperiodic orbits in a neighbourhood of the latter. The valleys appear when some orbits are locked inside a resonance. For example, looking at Fig. 10 for q=0.7 we observe a single peak very close to x₀=0.06 and, for , we see two peaks at both sides of a valley. In the first case the peak is due to the fact that we are crossing an island chain through the hyperbolic point, the orbits never fall inside the resonance. On the other hand, when we cross the island chain through the center of one island, i.e. through the elliptic point, we intersect twice the separatrix and some orbits are trapped by the resonance. The width of a peak or a valley is then a measure of the actual size of the resonance.

This can be visualized also from Fig. 11a where we plot a high-resolution surface of section for loops in the neighbourhood of the unstable periodic orbit at $x_0 \approx 0.06$ (4:9 resonance, see Sect. 6). The arrows in this figure indicate some of the resonances observed in the corresponding Fig. 10. In the latter figure, for we have marked two small peaks that should correspond to unstable periodic orbits but in Fig. 11a their presence is not evident. Nevertheless for $x_0\approx0.055$ one can distinguish a rather narrow resonance that should be responsible of one of the mentioned peaks (in fact, to that located more distant to the 4:9 resonance). A similar picture to that given by $\sigma_{\rm ls}$ comes from $\bar{{\cal J}}$ where Fig. 11b is representative to illustrate the MEGNO's behaviour. Here we plot the final value of $\bar{{\cal J}}$ for the same window shown in Fig. 11a and only for $\sim 400$ orbits with the same resolution in x₀ like that in Figs. 9 and 10. Along this interval $\bar{{\cal J}}$ is very close to 2 but we can appreciate the resonances observed in the surface of section as well as those marked in Fig. 10. In any case $1.98<\bar{{\cal J}}<2.035$ and the dotted line corresponds to the theoretical value for stable quasiperiodic motion, $\bar{\cal J}=2$.

Figure 10 as well as Fig. 12, where we plot the full domain of box family in separated windows (but with the same resolution in p_x0than in Fig. 9), are very illustrative to see how resonances in the neighbourhood of the stochastic layer overlap as the perturbation increases, leading to an enlargement of the layer and therefore to a larger domain of irregular motion. A few small and not too small islands are embedded in this chaotic zone, some of them are indicated in both figures. All the thin peaks observed in both figures correspond to hyperbolic orbits. Therefore we obtain a picture of the hyperbolic structure of the phase space that announces the future appearance of irregular motion as soon as we increase the perturbation (to be precise, irregular motion certainly exists around all these resonances but it is confined to a set of negligible measure). Additionally $\sigma_{\rm ls}$ reveals some details about the internal structure of the secondary resonances.