# The Fast Norm Vector Indicator (FNVI) method: A new dynamical parameter for detecting order and chaos in Hamiltonian systems

## Abstract

In the present article, we introduce and also deploy a new, simple, very fast and efficient method, the Fast Norm Vector Indicator (FNVI) in order to distinguish rapidly and with certainty between ordered and chaotic motion in Hamiltonian systems. This distinction is based on the different behavior of the FNVI for the two cases: the indicator after a very short transient period of fluctuation displays a nearly constant value for regular orbits, while it continues to fluctuate significantly for chaotic orbits. In order to quantify the results obtained by the FNVI method, we establish the dFNVI, which is the quantified numerical version of the FNVI. A thorough study of the method’s ability to achieve an early and clear detection of an orbit’s behavior is presented both in two and three degrees of freedom (2D and 3D) Hamiltonians. Exploiting the advantages of the dFNVI method, we demonstrate how one can rapidly identify even tiny regions of order or chaos in the phase space of Hamiltonian systems. The new method can also be applied in order to follow the time evolution of sticky orbits. A detailed comparison between the new FNVI method and some other well-known dynamical methods of chaos detection reveals the great efficiency and the leading role of this new dynamical indicator.

###### Keywords:

Hamiltonian systems; ordered and chaotic motion; new dynamical indicators^{1}

∎

## 1 Introduction

The issue of knowing whether the orbits of a dynamical system are ordered or chaotic is fundamental for the understanding of the behavior of the system in an extended area of modern science. In the dissipative case, this distinction is easily made as both types of motion are attracting. On the other hand, in conservative systems, distinguishing between regular and chaotic motion is often a very delicate and difficult issue (i.e. when the chaotic or ordered regions are small), especially in dynamical systems with many degrees of freedom, where one cannot easily visualize and interpret directly the behavior of the orbits. For this reason, it is often of great importance to possess fast and accurate tools in order to determine if an orbit is ordered or chaotic, independent of the dimension of the phase space of the dynamical system.

Let’s recall and present some well-known, classical methods that try to give an answer to the issue of determining the nature of an orbit.

(a) The inspection of the consequents of an orbit on a Poincaré surface of section (PSS). For 2D dynamical systems this technique has been used extensively, despite the problem of establishing a proper Poincaré surface of section in each case. However, the inspection becomes very difficult and also greatly deceiving in the case of dynamical systems with multidimensional phase space.

(b) The maximal Lyapunov Characteristic Exponent (LCE) of an orbit informs us whether an orbit is ordered or chaotic. If then the corresponding orbit is chaotic. Over thirty years ago, Benettin et al [2] studied theoretically the problem of the computation of all LCEs and proposed an algorithm for their numerical computation. In particular, is computed as the limit for of the quantity

(1) |

where and are the deviation vectors for a given orbit, at times and respectively. The time evolution of is given by solving the so-called variational equations. Generally, for almost all choices of initial deviations , the limit for of equation (1) gives always the same value of

(2) |

In practice, of course, since the exponential growth of occurs for short time intervals, one stops the evolution of after some time , records the computed , normalize the vector and repeats the calculation for another time interval , etc obtaining finally the as an average over many , as

(3) |

The basic drawback of the computation of is that, after every , the calculation starts from the beginning and may yield an altogether different than the interval. Thus, since is influenced by the whole evolution of , the time needed for (or the ) to converge is not known a priori and may become extremely long. This makes it often very difficult to determine whether finally tends to a positive value (chaos) or converges to zero (order). The main advantage of the LCE is that it can be applied easily to dynamical systems of any number of degrees of freedom.

(c) The frequency analysis method proposed by Laskar [14-16, 19, 20], which is based on the calculation of the basic frequencies of an orbit over a fixed interval of time. For orbits on the KAM tori, these frequencies are very accurate approximations of the actual frequencies of the dynamical system, but for chaotic orbits the computed values vary significantly in time and space. The frequency analysis method can be applied to systems with many degrees of freedom.

(d) The study of the dynamical spectra of orbits [6, 21, 28, 34, 35]. The distribution of the values of a given parameter along the orbit has been proved a very reliable and powerful tool for the exploration of the properties of motion in Hamiltonian systems of two and three degrees of freedom. The reader can find more interesting information about the dynamical spectra in [35] and also in the references of this article.

Over the last years, several methods have been introduced in order to characterize the nature of an orbit by studying the evolution of the deviation vectors, some of which are discussed in section 4. In the present paper we introduce and use a new, fast and easy to compute indicator: the Fast Norm Vector Indicator (FNVI). We focus our attention on the method of the FNVI, performing a systematic and thorough study of its behavior in the case of autonomous Hamiltonian systems with two (2D) and three (3D) degrees of freedom. The FNVI was found to fluctuate significantly for chaotic orbits, while it displays a nearly constant value for ordered orbits.

This article is organized as follows: in section 2 we provide the definition of the FNVI and we present results distinguishing between ordered and chaotic motion in two and three degrees of freedom (2D and 3D) Hamiltonians, comparing also the efficiency of the FNVI with the computation of the corresponding LCE. In the same section we establish a numerical criterion, that is the dFNVI, in order to quantify the results obtained by the FNVI method. In section 3 we demonstrate the ability of the new method to reveal the detailed structure of the dynamics in the phase space in both 2D and 3D dynamical systems. In the following section, we apply the new method to follow the evolution of a sticky orbit in the two-dimensional dynamical system. In section 5 we conduct a detailed comparison between the new FNVI/dFNVI method and with some other relatively modern well-known methods of chaos detection. Finally, in section 6 we summarize our results and we present a discussion and also the conclusions of the present research.

## 2 The definition of the FNVI method

The basic idea behind the FNVI method is the introduction of a simple and easily computed quantity that clearly identifies the ordered or chaotic nature of an orbit. The FNVI is defined as

(4) |

where is the time, denotes the Euclidean norm of the vector , while and are the norm vectors for a given orbit, at times and respectively. In practice, we stop the evolution of after some time , we record the computed and then we repeat the calculation for another time interval , etc obtaining finally the FNVI as a summation over many , . Here we must point out, that is the predefined time step of the numerical integration which remains constant during the entire predefined total integration time .

In order to apply the FNVI method, we shall investigate the nature of orbits in a dynamical system of perturbed harmonic oscillators given by the 3D potential

(5) |

where is the common frequency of the oscillations along the , and axis, while is the strength of the perturbation. Potential (5) represents three coupled harmonic oscillators in the case of the 1:1:1 resonance. Potentials of this type are also known as perturbed elliptic oscillators [1, 3, 7]. The basic reason for the choice of potential (5) is that perturbed elliptic oscillators appear very often in galactic dynamics and also in atomic-particle physics [7]. A second reason is that is displays exact periodic orbits, interesting sticky orbits together with large unified chaotic domains. Therefore, it gives us a great opportunity to test and prove the effectiveness and the reliability of the new FNVI method.

The Hamiltonian function to the potential (5) reads

(6) | |||||

where , and are the momenta per unit mass conjugate to , and respectively, while is the numerical of the Hamiltonian.

The corresponding Hamiltonian for the 2D system can easily be obtained if we set in equation (6). Then

(7) | |||||

where is the numerical value of the 2D Hamiltonian.

The outcomes of the present research are mainly based on the numerical integration of the equations of motion

(8) |

where the dot indicates derivative with respect to the time. We integrated numerically the equations of motion with a double precision Bulirsch-Stöer algorithm in FORTRAN 95. The accuracy of our results was checked by the constancy of the energy integrals (6) and (7), which were conserved up to the fifteenth significant decimal point.

We keep the involved parameters of the two systems fixed at the values and , while the value of the energy or is treated as a parameter.

A simple qualitative way of studying the dynamics of a Hamiltonian system is by plotting the successive intersections of the orbits with a Poincaré surface of section [17]. This method has been extensively applied to 2D Hamiltonians, as in these systems the PSS is a two-dimensional plane. In 3D systems, however, the PSS is four dimensional and the behavior of the orbits cannot be easily visualized. One way to overcome this problem is to project the PSS to spaces with lower dimensions (see, [31, 32]). However, even these projections are often very complicated and difficult to be interpreted.

In order to illustrate the behavior of the FNVI in 2D and 3D dynamical systems, we first consider some representative ordered and chaotic orbits. In Figure 1a we plot the intersection points of an ordered and a chaotic orbit of the dynamical system (7), with a PSS defined by and . The points of the ordered orbit lie on a torus and form a smooth closed curve on the PSS. On the other hand, the points of the chaotic orbit appear randomly scattered. The outermost black solid line shown in Fig. 1a is the Zero Velocity Curve (ZVC). The equation of the limiting curve ZVC (that it the curve containing all the invariant curves of the 2D system) is defined by the equation

(9) |

The time evolution of the FNVI for these two orbits and for a time period of time units is plotted in Figure 1b. In the case of the ordered orbit the FNVI remains almost constant, while in the case of the chaotic orbit it displays large fluctuations. The initial conditions for the regular orbit are: , and , while the initial value of is found from the energy integral (7). The chaotic orbit has initial conditions: , and . For both orbits the value of the energy is , while the the PSS shown in Fig. 1a was integrated for a time period of 5 time units. In order to have a better and more closer view of the evolution of the FNVI for these two orbits, we present separately in Figure 2a-b the two diagrams. We observe, that in Fig. 2a which corresponds to the regular orbit, the FNVI after a small transient period of fluctuation for about 200 time units, it stabilizes at a value and remains almost constant. On the contrary, in the case of the chaotic orbit shown in Fig. 2b the FNVI displays a completely different character with large and abrupt fluctuations throughout the entire period of the numerical integration.

One may reasonably assume that the shape of the FNVI depends on the particular value of the total time of the numerical integration. Thus, if we integrate the FNVI of the chaotic orbit shown in Fig. 2b for a larger time period it might displays a similar pattern as the one of the ordered orbit shown in Fig. 2a. In order to clarify this issue, we computed the FNVI for the same two 2D orbits discussed in Fig. 2a-b but for a time period of time units. The results are presented in Figs. 3b and 3d respectively and compared with the corresponding LCE in order to test their validity. In Figure 3a we see the evolution of the LCE of the ordered orbit. As expected the value of the LCE tends to zero as the time evolves. The evolution of the FNVI is presented in Figure 3b. Again, the FNVI after a small transient period of fluctuation for about 200 time units, it stabilizes at a value and remains almost constant. This indicates the regular character of the orbit. On the other hand, the computation of the maximal LCE, using Eqs. (1) and (2), despite its usefulness in many cases, does not have the same convergence properties over the same time interval. This becomes evident in Figure 3c where we plot the evolution of the LCE for the chaotic orbit. The computation of the LCE up to or even up to time units, still shows no clear evidence of convergence. Although Fig. 3c suggests that the orbit might probably be chaotic since LCE remains different from zero for large time intervals, it does not allow us to conclude its chaotic nature with certainty, and so further computation of the LCE is needed. Thus, it becomes evident that an advantage of the FNVI, with respect to the computation of the LCE, is that the current value of the FNVI is sufficient to determine the chaotic nature of an orbit, in contrast to the maximal LCE, where the whole evolution of the deviation vector affects the computed value of LCE. Figure 3d depicts the evolution of the FNVI for the chaotic orbit. Therefore, it becomes clear that the shape of the FNVI does not depend on the time of the numerical integration. We observe in Fig. 3d that the FNVI displays once more large and abrupt fluctuations throughout the time evolution. Note that in Figs. 3a-d the horizontal axis corresponding to the time is on log scale, in order to visualize more clearly the evolution of the indicators throughout the entire time interval. As the criterion by which we judge using the FNVI method, whether an orbit is regular or chaotic is qualitative rather than quantitative, we could be sure for the chaotic nature of the orbit after only time units of numerical integration.

Let us now present some results about the 3D Hamiltonian system (6). In Figure 4a we observe the evolution of the LCE of an ordered 3D orbit with initial conditions: , while the initial value of is found from the energy integral (6). As expected, the value of the LCE tends to zero as the time evolves. The evolution of the FNVI is presented in Figure 4b. The FNVI after a small transient period of fluctuation for about 250 time units, it stabilizes at a value and remains almost constant, indicating the regular character of the orbit. In Figure 4c we see the evolution of the LCE of a chaotic 3D orbit. In this case the initial conditions of the chaotic orbit are: . The LCE remains different from zero, for a time period of time units, which implies the chaotic nature of the orbit. The corresponding evolution of the FNVI is presented in Figure 4d. It is clear that the pattern presented in Fig. 4d is completely different from that shown in Fig. 4b, where the 3D orbit is regular. For the 3D chaotic orbit, the FNVI displays large fluctuations and a highly irregular shape. Note that the axis in Figs. 4a-d is in log scale. For both 3D orbits the value of the energy is .

Here, we have to point out that the criterion by which we determine so far the regular or chaotic nature of an orbit (2D or 3D), using the FNVI method is purely qualitative. In the case of a regular orbit the FNVI after a small transient period of fluctuation it stabilizes at a value and remains almost constant. On the other hand, when the orbit is chaotic the evolution of the corresponding FNVI follows a completely different pattern. It does not shows any evidence of convergence and displays large and abrupt fluctuations. This criterion has been obtained by testing a large number of regular and chaotic orbits in both dynamical systems (2D and 3D). Of course, since our criterion is qualitative we have to inspect the shape of FNVI each time in order to characterize an orbit. Obviously, this is not very practical when someone wants to check a large volume of orbits, so as to form an idea about the global structure of the dynamical system. Thus, we need to establish a new numerical criterion, in order to quantify the results obtained by the FNVI method. This criterion can be derived by looking the shape of the FNVI shown in Figs. 3b, 3d, 4b and 4d. One may observe, that when the orbit is regular the FNVI remains almost constant, while in the case of a chaotic orbit it displays high fluctuations. We are going to exploit this significant difference in order to obtain our quantitative criterion. Thus, we calculate the maximum and the minimum value of FNVI when . We take this particular time interval because in the case where the orbit is regular, for time units there is a transient period of fluctuation, which may create a malfunction to our criterion. Using the above procedure we can define the

(10) |

where and are the maximum and the minimum value of FNVI respectively when . The value of dFNVI can provide us the quantitative criterion that we seek. We note, that it is not easy to define a threshold value, so that the dFNVI being larger than this value reliably signifies chaoticity. Nevertheless, extensive numerical experiments in both dynamical systems indicate that in general, a good guess for this value could be 0.05. Thus, when dFNVI 0.05 the orbit is chaotic, while when dFNVI 0.05 the orbit ir ordered. This threshold value applies in both 2D and 3D orbits. For the 2D regular orbit presented in Fig. 2a dFNVI = 0.009, while for the 2D chaotic orbit of Fig. 2b dFNVI = 0.45. Similarly, for the 3D regular orbit presented in Fig. 4b dFNVI = 0.006, while for the 3D chaotic orbit of Fig. 4d dFNVI = 0.52. Therefore, we may conclude that the quantification of the FNVI method has been archived successfully. In the next section, we shall provide more extensive and detailed results using the dFNVI.

## 3 Distinguishing between regions of order and chaos

The dFNVI offers indeed a very easy and efficient method for distinguishing the chaotic versus ordered nature of orbits in a variety of problems. In the present section, we use it for identifying regions of phase space where large scale ordered and chaotic motion are both present. In particular, we shall apply the dFNVI method in order to reveal the structure in both Hamiltonian systems of two (2D) and three (3D) degrees of freedom.

### 3.1 Order and chaos in the 2D dynamical system

In Figure 5a-d we present detailed plots of the , , PSS of the 2D dynamical system (7) for different values of the energy . Fig. 5a shows the PSS plot, when . One can see, that regions of ordered motion around stable periodic orbits are seen to coexist with chaotic regions filled by scattered points. In particular, a large part of the phase plane is covered by chaotic orbits, while the regular regions are occupied by invariant curves mainly around the points and . The above points give the position of two exact periodic orbits. The first is the axis, while the second describes the straight-line orbits on the phase plane. Furthermore, one observes a large number of smaller islands of invariant curves produced by secondary resonances and some sticky regions as well. Fig. 5b is similar to Fig. 5a but when . Here, we see that the majority of the phase plane is covered by a unified chaotic sea. The regular regions are confined only around the straight line periodic orbits, while the axis has now become unstable. Some smaller islands of invariant curves are also present. The most interesting feature observed in this case, is the sticky regions around the chain of islands of invariant curves on the as well as on the axis. If we increase the value of energy to , we obtain the phase plane shown in Fig. 5c. In this case, the axis has returned to stability, while the straight line periodic orbits are unstable. Almost all the phase plane is chaotic, except of a small region around the center and some small islands of invariant curves which are products of secondary resonant orbits. A sticky region around the chain of the small islands of invariant curves near the center is also observed. Fig. 5d is similar to Fig. 5a, but when . Here, almost the entire phase plane is covered by chaotic orbits. There is a very small regular region confined near the center. A careful observation also reveals some very tiny regular islands of invariant curves embedded in the vast chaotic sea.

In order to demonstrate and prove the effectiveness of the dFNVI method, we first consider orbits whose initial conditions lie on the lines and . In particular, we take 5 equally spaced initial conditions on these lines and we compute the value of the dFNVI for each one. The results are presented in Figure 6a-h where we plot the dFNVI as a function of the and coordinates of the initial conditions of these orbits for time units. In all panels the data points are line connected, so that the changes of the dFNVI values are clearly visible. With green dots we represent the initial conditions or corresponding to regular orbits, while the initial conditions producing chaotic orbits are marked with red dots. Note that there are intervals where the dFNVI has low values (e.g. lower than 0.05), which correspond to ordered motion inside the islands of stability crossed by the and lines shown in Fig. 5a-d. There also exist regions where the dFNVI has large values (e.g. larger than 0.1) denoting that in these regions the motion is chaotic. These intervals correspond to the regions of scattered points crossed by the and lines in Fig. 5a-d. Although most of the initial conditions give large or very small values for the dFNVI, there also exist initial conditions that have intermediate values of the dFNVI (0.05 dFNVI 0.1). These initial conditions correspond to sticky chaotic orbits, remaining for long time intervals at the borders of islands, whose chaotic nature will be revealed later on. In Figs. 6a and 6b we observe the values of the dFNVI computed for time units for orbits with initial conditions on the and lines on the PSS shown in Fig. 5a, as a function of the coordinate and the momenta respectively. Similarly, Figs. 6c and 6d correspond to the PSS of Fig. 5b, Figs. 6e and 6f correspond to the PSS of Fig. 5c and Figs. 6g and 6h correspond to the PSS shown in Fig. 5d.

In Fig. 6a, around there exists a group of points inside a large chaotic region having dFNVI 0.05. These points correspond to orbits with initial conditions inside a small stability island, which is not very visible in the detailed PSS plot of Fig. 5a. Also the points with 0.91 in Fig. 6b, have been characterized by the dFNVI as initial conditions corresponding to chaotic orbits, while all their neighboring points have dFNVI 0.05. These points actually correspond to a weak chaotic separatrix inside the vast domain of stability, which can be revealed only after a very high magnification of this region of the PSS. So, we see that the systematic application of the dFNVI method can reveal very fine details of the structure in a dynamical system.

It will be of particular interest, to calculate the mean value of dFNVI for the regular and chaotic orbits discussed in Fig. 6a-h. In other words, we are going to calculate the mean value of dFNVI for the regular orbits, that is and the mean value of dFNVI for the chaotic orbits , with initial conditions along the and the axis for different values of the energy . In Table 1 we present our results. It is evident, that as the value of the energy increases, we have an increase at the mean value of the dFNVI for both regular and chaotic orbits. In the case of chaotic orbits, the increase of indicates that the chaoticity of the dynamical system increases as the value of the energy is amplified. We shall come back to this point later in this section.

- axis | - axis \bigstrut | |||
---|---|---|---|---|

1.0 | 0.00725 | 0.32251 | 0.00551 | 0.28687 |

2.0 | 0.01136 | 1.14872 | 0.01922 | 1.15079 |

4.0 | 0.01732 | 1.27763 | 0.02317 | 1.21792 |

7.2 | 0.02181 | 1.43602 | 0.02761 | 1.42706 |

By carrying out the previously presented analysis for initial conditions not only along a line but on the whole phase plane of a PSS and giving to each point a color according to the value of the dFNVI, we can have a clear picture of the regions where chaotic or ordered motion occurs. The outcomes of this procedure for the 2D dynamical system (7), using a dense grid of initial conditions on the PSS, are presented in Figure 7a-d. The value of the energy in Figs. 7a, 7b, 7c and 7d is the same as in Figs. 5a, 5b, 5c and 5d respectively. Thus, in Figs. 7a-d we clearly distinguish between green regions, where the motion is ordered and red regions, where it is chaotic. The outermost black solid line shown in Figs. 7a-d is the limiting curve (ZVC) defined by Eq. (9). It is worth mentioning, that in Fig. 7a we can observe small islands of stability inside the large chaotic sea, which are not very visible in the detailed PSS of Fig. 5a, such as that for 1.05, 0. Although all the orbits with initial conditions in Figs. 7a-d were computed for only = time units (such as Figs. 6a-h), this time was sufficient enough for the clear revelation of small ordered regions inside the chaotic domains. We must point out, that the outcomes derived using the dFNVI method regarding the structure of the phase plane coincide with those obtained using the PSS technique. In other words, we see that the structure of the phase planes shown in Figs. 7a-d are practically identical with those presented in Figs. 5a-d. For a grid of about equally spaced initial conditions , we need about 7 h of CPU time on a Pentium Dual Core Processor at 2.2GHz PC, in order to construct each grid-plot shown in Figs. 7a-d.

Grid-plots such as those of Fig. 7a-d, apart from presenting the regions of order and chaos, can also be used to provide very accurate estimations regarding the fraction of the phase-space volume occupied by chaotic or ordered orbits and also give us good initial guesses for the location of stable periodic orbits, in regions where the motion is ordered. One can conclude form Figs. 5a-d or from Figs. 7a-d, that the fraction of the phase plane covered by chaotic orbits increases as the value of the energy increases. Figure 8 shows a plot of the percentage of the phase plane covered by chaotic orbits versus . We observe that the percentage increases rapidly, as the value of increases. Dots indicate the values of the percentage obtained numerically from the grid-plots, while the solid line is a fourth degree polynomial fitting curve. We must point out that the percentage is calculated as follows: when constructing the grid-plots, we count the initial conditions corresponding to regular orbits and also the initial conditions corresponding to chaotic orbits. Then, we divide the number of chaotic orbits by the total number of orbits and thus, we obtain the percentage for each value of the energy . Note, that in every case the total number of the calculated orbits is about , while the exact number of chaotic orbits varies depending on the value of the energy .

In order to have an estimation of the degree of chaos of the dynamical system from another point of view, we have plotted the average value of the maximal LCE versus . The results are shown in Figure 9a. Note, that the LCE increases linearly with . Here, we must point out that it is well known that the value of the LCE is different in each chaotic component [22]. As we have in all cases regular regions and only one unified chaotic sea in each phase plane, we calculate the average value of LCE by taking orbits with different and random initial conditions in the chaotic domain for every value of the energy . Note that all calculated LCEs are different on the fifth decimal point in the same chaotic region. We follow the same procedure using now the dFNVI. For the same initial conditions corresponding to chaotic orbits for each value of the energy , we calculated the dFNVI . The results are presented in Figure 9b. Again the dFNVI increases linearly with . The increase of dFNVI indicates that the chaoticity of the dynamical system increases as the value of the energy is amplified. As the behavior of dFNVI coincides with the behavior of LCE , we prove once more the reliability of dFNVI method. However, the main advantage of the dFNVI method is that it requires only integration time units, while the LCE needs at least about 5 integration time units so as to provide reliable and definitive evidence. Thus, the dFNVI is a much more faster indicator than the LCE.

### 3.2 Order and chaos in the 3D dynamical system

In the case of 3D Hamiltonian systems the PSS is now four dimensional and thus, not so useful as in the 2D systems. On the other hand, the dFNVI method has the ability once more to identify successfully regions of order and chaos in the phase space. To prove this, let us start with initial conditions on a 4D grid of the PSS. In this way, we find again regions of order and chaos, which may be visualized, if we restrict our study to a subspace of the whole 6D phase space. We consider orbits with initial conditions , , while the initial value of is always obtained from the energy integral (6). In particular, we define a value of which is kept constant and then we calculate the dFNVI of 3D orbits with initial conditions , . Thus, we are able to construct again a 2D plot depicting the plane but with an additional value of . All the initial conditions of the 3D orbits lie inside the limiting curve defined by

(11) |

Our results are presented in Figure 10 a-d, where we give four grid-plots for the same value of the energy , but for different initial value of . Again, we can distinguish between regions of ordered (colored in green, where dFNVI 0.05) and chaotic motion (colored in red, where dFNVI 0.05). In Fig. 10a where we see that the structure of the plane is quite similar with the corresponding 2D grid-plot shown in Fig. 7a. We also see that well defined islands of stability inside the unified chaotic sea still exist. The main difference from the 2D grid-plot of Fig. 7a is that now in the 3D system, a large portion of orbits around the stable periodic points on the axis have altered their nature form ordered to chaotic. Fig. 10b is similar to Fig. 10a but when . In this case, the structure of the plane is very similar to that shown in Fig. 10a. One may observe, that the percentage corresponds to chaotic orbits has been increased and moreover the islands of stability have begun to destabilize and lose their well defined structure. In Fig. 10c we increase the initial value of to 0.3. In this case, it is evident that the initial conditions correspond to ordered orbits are delocalized, since now we can see no signs of well defined islands of stability. Finally, in Fig. 10d we present a grid-plot when . Here, the vast majority of the initial conditions correspond to chaotic orbits, while the initial conditions correspond to ordered orbits are completely delocalized and randomly scattered mainly in the central part of the plane. All the 3D orbits with initial conditions in Figs. 10a-d were computed for only time units. With a closer look to Figs. 10a-d one may conclude the following two main points: (i) the increase of the initial value of has as a result the increase of the percentage on the plane corresponding to chaotic orbits and (ii) as we amplify the value of we observe that the entire area of the plane defined by Eq. (11) is reduced and becomes smaller. This can easily be explained. The coordinates and the momenta of the 3D orbits are connected through the energy integral (6) and since we increase the value of while keeping constant the value of the energy this has as the result the permissible values of the initial conditions to be reduced. For a grid of about equally spaced initial conditions , we need about 7.5 h of CPU time on a Pentium Dual Core Processor at 2.2GHz PC, in order to construct each grid-plot shown in Figs. 10a-d.

The grid-plots shown in Figs. 10a-d could be considered as “phase planes” but not with the strict definition. In fact, these plots show the structure of the subspace for a given value of the . Therefore, we study the subspace , , , with of the 6D phase space of the 3D dynamical system. Note that grid-plots such as those of Figs. 10a-d, apart from presenting the initial conditions correspond to ordered and chaotic orbits, can also be used to estimate the fraction of the phase-space volume occupied by chaotic or ordered 3D orbits. In Figure 11, we present a plot depicting the relationship between the percentage of initial conditions corresponding to chaotic 3D orbits and the initial value of . We observe that the chaotic percentage increases linearly with . Our numerical calculation suggest, that when for eventually all the initial conditions correspond to chaotic orbits, while regular motion, if any, is negligible. Note, that the plot shown in Fig. 11 corresponds to . For different values of the energy the relationship between and is quite similar, but as we amplify the value of the energy , the increase of the chaotic percentage is more rapid and does not follow a linear pattern any longer.

Finally, in order to have an estimation of the degree of chaos of the 3D dynamical system from another point of view, we have plotted the average value of the maximal LCE versus the value of the energy . The results are shown in Figure 12a. Note that the LCE increases linearly with , when the initial value of is 0.01. As we have regular regions and only one unified chaotic sea in the planes (see Fig. 10a), we calculate the average value of the LCE by taking orbits with different and random initial conditions , and in the chaotic domain for every value of the energy . Note, that all calculated LCEs are different on the fifth decimal point in the same chaotic region. For different initial value of the , our numerical experiments indicate that the relationship between LCE and remains linear and similar to that shown in Fig. 12a. Here we should point out, that the LCE of the 3D system is smaller than the corresponding LCE of the 2D system shown in Fig. 9a. We follow the same procedure using now the dFNVI method. For the same initial conditions , and corresponding to chaotic 3D orbits for each value of the energy , we calculate the dFNVI . The results are presented in Figure 12b. Again, the dFNVI increases linearly with . The increase of the dFNVI indicates that the chaoticity of the dynamical system increases as the value of the energy is amplified. For different initial value of the , our numerical calculations indicate that the relationship between dFNVI and also remains linear and similar to that shown in Fig. 12b. As the behavior of dFNVI coincides with the behavior of LCE , we prove the reliability of dFNVI method also in the 3D dynamical system. Of particular interest is the fact, that the dFNVI of the 3D system is also smaller than the corresponding dFNVI of the 2D system shown in Fig. 9b. This strongly suggests that the chaoticity in the dynamical system of three degrees of freedom(3D) is smaller that in the corresponding dynamical system of two degrees of freedom (2D), when .

## 4 Application in the case of a sticky orbit

In this section, we shall try to test the FNVI method’s ability posing a somehow more difficult task, that is to distinguish between ordered and chaotic motion, in the case of chaotic orbits that seem to be regular for a number of periods. Such orbits, called “sticky orbits” and they usually grow near the outer regions of an island of stability [18]. Cantori that exist outside that region may restrict the motion in a very thin layer for a long time interval. After that time interval, the orbit can escape to the surrounding chaotic region, through the gaps of the cantorus. Due to this peculiar behavior, a sticky orbit often resembles an ordered one seen in a Poincaré surface of section, until the orbit escapes to the chaotic domain.

We follow the evolution of a two-dimensional sticky orbit which produces two sets of five small islands of invariant curves near the axis, using the PSS technique. This sticky orbit has initial conditions: , and , while the value of is obtained from the energy integral (2). Our results are presented in Fig. 13 (a-d). Fig. 13a shows the regions formed in the phase plane by this orbit, for a period of 1200 time units of numerical integration. We see two sets of five separate islands of invariant curves. The outermost black solid line shown in Fig. 13a is the Zero Velocity Curve (ZVC) defined by Eq. (9). When time units (Fig. 13a), one can see that the orbit is chaotic merely zooming in the plot and seeing that the points are not smoothly distributed over a well defined one-dimensional curve. In Fig. 13b the same orbit was calculated for 1700 time units. Now, we observe that in each case the five small islands of invariant curves were joined together. By observing Figs. 13a and 13b, it is obvious that from the pictures presented in the PSS no safe conclusion can be reached concerning the regular or chaotic behavior of this orbit. Fig. 13c depicts the same orbit calculated for a time period of 12500 time units. It is evident, that now the test particle has left the sticky region and begins to enter the vast surrounding chaotic region. In Fig. 13d we have calculated the orbit for 5 time units. The chaotic nature of the orbit is now completely revealed, since the scattered points in the phase plane fill the entire available chaotic domain.

A better and more enlightening view for the evolution of the sticky orbit can be seen using the dynamical spectrum [33,35]. Fig. 14a shows the spectrum of the sticky orbit for 1150 time units of numerical integration. Here we observe ten separate -type structure with a large number of additional small and large peaks indicating the sticky nature of the orbit. In Fig. 14b we present the spectrum computed for a time interval of 1670 time units. In this case, the ten structures have joined together producing two separate new structures. Fig. 14c depicts the spectrum of the same orbit calculated for 12200 time units of numerical integration. One may observe, that the two separate structures shown in Fig. 14b have now joined together producing a single spectrum. The shape of this spectrum has all the characteristics of a chaotic spectrum and thus strongly indicates that the test particle has undoubtedly left from the sticky region and continued its journey to the surrounding chaotic sea. Fig. 14d shows the spectrum computed for 5 time units. We see that the spectrum has the usual shape of a spectrum corresponding to a chaotic orbit. Here we must point out, that the time intervals regarding the gradual time evolution of the sticky orbit obtained by the spectrum are very close to those derived by the formation of the regions in the phase plane, shown in Fig. 13 (a-d).

Now let’s use some other dynamical indicators in order to follow the evolution of the sticky orbit and also compare and verify our results. Fig. 15 (a-b) shows the time evolution of the LCE and the FNVI for the same sticky orbit discussed and studied earlier. Fig. 15a shows the time evolution of the LCE. It is evident that after 12600 time units the orbit gradually becomes chaotic. In Fig. 15b we present the time evolution of the FNVI for this orbit. We observe, that the pattern of the time evolution of the FNVI is almost identical with the pattern of the LCE shown in Fig. 15a. Furthermore, the FNVI clearly indicates that after 12100 time units the orbit gradually alters its nature to chaotic. Note that in Fig. 15a-b the axis is in log scale in order to have a more clear view of the time evolution. Therefore, we may conclude that the results obtained by the PSS technique, the dynamical spectrum, the LCE and also the FNVI method shown in Figs. 13a-d, 14a-d, 15a and 15b respectively, coincide that the sticky orbit becomes chaotic after a time interval of about 12500 time units of numerical integration. Thus, we may conclude that the new FNVI method has the ability to follow the time evolution of a two-dimensional sticky orbit and provide reliable results regarding the time needed for the transition from sticky motion to chaos.

## 5 Comparison with other dynamical indicators

The results presented in the previous sections reveal, that the dFNVI is a simple, efficient and easy to compute tool for the distinguishing between ordered and chaotic motion in Hamiltonian systems. Implementing the dFNVI method is a very easy computational task, as we only have to follow the evolution of the norm of an orbits’s vector throughout the integration time interval and compute the maximum and the minimum value of FNVI when . In the case of regular motion the FNVI after a small transient period of fluctuation it settles down to a value and remains almost constant, while the dFNVI is always smaller than 0.05. In particular, our extensive numerical calculations suggest, that when an orbits is regular the dFNVI is much smaller that 0.01. However, we decided to increase the threshold value to 0.05 in order to obtain more secure and reliable results. On the other hand, in the case of chaotic motion the FNVI continues to fluctuate randomly without giving any sign of convergence, while the numerical criterion, that is the dFNVI is always larger than 0.05. It is exactly this different behavior of the FNVI regarding the convergence, that makes it an ideal tool of chaos detection. However, the results obtained by the FNVI are only qualitative. Thus, we established a numerical criterion the dFNVI, which can provide quantitative results about the nature of an orbit. As we have seen in the previous sections, the results obtained by the dFNVI are completely reliable, conclusive and beyond any doubt. The dFNVI helps us decide the chaotic or regular nature of orbits faster and with less computational effort than the estimation of the maximal LCE. This happens because the time needed for the LCE in order to give a clear and undoubted indication of convergence to non-zero values is usually much greater than the time needed for the dFNVI to provide reliable evidence regarding the character of an orbit.

Many other dynamical indicators have been introduced in the recent years, some of which are compared in this section with the dFNVI method. In order to verify and prove, once more, the effectiveness and the reliability of our new method, we shall compare our results with four other well known dynamical indicators: (a) The Fast Lyapunov Indicator (FLI), (b) The Relative Lyapunov Indicator (RLI), (c) The Smaller Alignment Index (SALI) and (d) The Generalized Alignment Index (GALI). At this point we believe, that it should be wise to recall the definitions of these indicators, in order to help the reader to compare the results of each method.

The FLI introduced in [8, 9] has been proved a very fast, reliable and effective tool, which can be defined as

(12) |

where is a deviation vector. The computation of the FLI on a relatively short time is enough to discriminate between chaotic and regular orbits. The FLI of a regular orbit increases linearly, while for a chaotic orbit, the FLI increases super-exponentially. Moreover, the FLI may be used to discriminate among regular motion between non resonant and resonant orbits. In order to avoid the linear trend of the FLI along regular orbits, we can use the de-trended FLI (DFLI), which is simply the FLI divided by . Thus,

(13) |

In this case, the classification of an orbit is based on the clearly distinct behavior of the DFLI pattern. In particular, in the case of a regular orbit the DFLI is almost constant with values lower than 10, while in the case of a chaotic orbit, we observe a very rapid increase of the DFLI. We believe, that using definition (13) instead of (12) we can obtain better and more conclusive results. The DFLI is a relatively fast indicator, as it needs only about 5 times units in order to provide reliable results regarding the nature of an orbit. Note, that for the calculation of the DFLI, we need the assistance of the variational equations and of course the computation of the deviation vector .

Another interesting dynamical indicator is the RLI proposed by Sándor et al [23]. The RLI is practically the absolute value of the difference of of two initially nearby orbits and can be defined as

(14) |

where

(15) |

while . Very small values of the RLI denote ordered motion, while larger values denote chaotic behavior (for more information on the RLI method see [23]). The RLI method needs at least about 5 times units in order to provide reliable results regarding the nature of an orbit. We must point out, that the computation of the RLI requires the time evolution of two orbits and two deviation vectors ( and one for each orbit), while the computation of the dFNVI is much faster as we compute only the evolution of the main orbit without the use of any deviation vector.

Using the variational equations and the evolution of the deviation vectors we can compute the SALI [25, 26]. The basic idea behind the SALI method is the introduction of a simple quantity that clearly indicates if a deviation vector is aligned with the direction of the eigenvector which corresponds to the maximal LCE. In general, any two randomly chosen initial deviation vectors and will become aligned with the most unstable direction and the angle between them will rapidly tend to zero [2]. Thus, we check if the two vectors have the same direction in phase space, which is equivalent to the computation of the above-mentioned angle. More specifically, we follow simultaneously the time evolution of an orbit with initial condition and two deviation vectors with initial conditions and . As we are only interested in the directions of these two vectors we normalize them, at every time step, keeping their norm equal to 1. This controls the exponential increase of the norm of the vectors and avoids overflow problems. Since, in the case of chaotic orbits the normalized vectors point to the same direction and become equal or opposite in sign, the minimum of the norms of their sum (antiparallel alignment index) or difference (parallel alignment index) tends to zero. So, the SALI is defined as

(16) |

where and are the parallel and the antiparallel alignment indices respectively. From the above definition it is evident that SALI(t) and when SALI = 0 the two normalized vectors have the same direction, being equal or opposite. Implementing the SALI method is an easy computational task, as we only have to follow the evolution of an orbit and of two deviation vectors, computing in every time step the minimum norm of the difference and the addition of these vectors. In the case of chaotic motion the SALI eventually tends exponentially to zero, reaching rapidly very small values or even the limit of the accuracy of the computer . On the other hand, in the case of ordered motion the SALI fluctuates around non-zero values. The SALI has a clear physical meaning as zero, or a very small value of the index, signifies the alignment of the two deviation vectors. An advantage of the method is that the index ranges in a defined interval (SALI and so very small values of the SALI (e.g., smaller than ) establish the chaotic nature of an orbit beyond any doubt. The SALI is a very fast indicator, as it needs only about times units in order to reveal the regular or chaotic character of an orbit.

Recently, Skokos et al [27] generalized and improved considerably the SALI method mentioned above, by introducing the GALI method. This indicator retains the advantages of the SALI i.e. its simplicity and efficiency in distinguishing between regular and chaotic motion but, in addition, is faster than the SALI, displays power law decays that depend on torus dimensionality and can also be applied successfully to cases where the SALI is inconclusive, like in the case of chaotic orbits whose two largest Lyapunov exponents are equal or almost equal. For the computation of the GALI we use information from the evolution of more than two deviation vectors with respect to the reference orbit, while SALI’s computation requires the evolution of only two such vectors. In particular, is proportional to “volume” elements formed by initially linearly independent unit deviation vectors whose magnitude is normalized to unity at every time step. If the orbit is chaotic, goes to zero exponentially fast by the law

(17) |

If, on the other hand, the orbit lies in an -dimensional torus, displays the following behaviors: Either

(18) |

or, if , it decays with different power laws, depending on the number of deviation vectors which initially lie in the tangent space of the torus, i.e.,

(19) |

So, the GALI method allows us to study more efficiently the geometrical properties of the dynamics in the neighborhood of an orbit, especially in higher dimensions, where it allows for a much faster determination of its chaotic nature, overcoming the limitations of the SALI method. In the case of regular motion, is either a constant, or decays by power laws that depend on the dimensionality of the subspace in which the orbit lies, which can prove useful e.g., if our orbits are in a “sticky” region, or if our system happens to possess fewer or more than independent integrals of the motion (i.e. is partially integrable or super-integrable respectively). It should be mentioned, that both the SALI and the GALI methods are very reliable and efficient. However, in order to apply them we have to use the variational equations and also several sets of deviation vectors. On the contrary, the dFNVI method is much faster than these methods because it does nor require the use of the variational equations and the deviation vectors. Simply we follow the time evolution of an orbit with initial condition and computing its norm by integrating only the basic equations of motion.

In what follows, we shall present and compare representative examples for 2D and 3D orbits of the dynamical systems (7) and (6) respectively, using the above four indicators. In Figure 16a-d, we can see the results provided by the four indicators for a regular 2D orbit with initial conditions: , and , while the initial value of is found from the energy integral (7). The value of the energy is . In Fig. 16a we present the evolution of the DFLI for a time period of 5 time units. We observe that the DFLI fluctuates around small values (DFLI 10), which indicates regular motion. The initial deviation vector, , used for the orbit shown in Fig. 16a is . Fig. 16b depicts the evolution of the RLI for a time interval of time units, for the same regular orbit. The RLI clearly reveals the ordered nature of the orbit, since RLI throughout the time evolution. The initial deviation vectors, used for the computation of the RLI are and , while . Note that the axis is in log scale. In Fig. 16c we present the evolution of SALI for a time interval of 5 time units. Also this method indicates the regular character of the orbit, since it exhibits small fluctuations around non-zero values. The initial deviation vectors used for the orbit of Fig. 16c are and , but in general, any other choice of initial values leads to similar behavior of the SALI. From Eqs. (18) and (19) it follows that in the case of a Hamiltonian system with = 2 degrees of freedom will always remains different from zero, while and should decay to zero following a power law, whose exponent depends on the number of the deviation vectors that are initially tangent to the torus on which the orbit lies. Now, for the regular orbit of the 2D Hamiltonian (7) and for a random choice of initial deviation vectors, we expect the indices to behave as

(20) |

Fig. 16d shows the behavior of for the same regular 2D orbit. Indeed, the indices obey to the power laws respectively given in Eq. (20). In Fig. 16d with red color we plot the power law , while the power law is colored in green. Note that the axis is in log scale. The evolution of the FNVI for this regular 2D orbit, for a time interval of time units, is presented in Fig. 3a, while dFNVI = 0.009. We see, that all the outcomes derived using five different dynamical indicators coincide to the ordered nature of the orbit.

Figure 17a-d is similar to Fig. 16a-d but for a chaotic 2D orbit. In Fig. 17a-d, we present the results provided by the four dynamical indicators for a chaotic 2D orbit with initial conditions: , and , while the initial value of is found from the energy integral (7). The value of the energy is . In Fig. 17a we present the evolution of the DFLI for a time period of 5 time units. We observe that the DFLI displays an exponential growth with time, which indicates chaotic motion. The initial deviation vector used for the orbit shown in Fig. 17a is . Fig. 17b depicts the evolution of the RLI for a time interval of time units, for the same chaotic orbit. Approximately for the first 350 time units, the RLI has low values (RLI ) but for the rest time interval RLI . Note that the axis is in log scale. Therefore, we may conclude that the orbit is chaotic. Note, that the RLI method gives inconclusive or even misleading results for very small time intervals of the numerical integration. The initial deviation vectors, used for the computation of the RLI are and , while . In Fig. 17c we present the evolution of the SALI. We see that after a small transient time, the SALI falls abruptly to zero. At time units the SALI becomes zero, as it has reached the limit of the accuracy of the computer , which means that the two deviation vectors have the same direction. Thus, after time units the two normalized vectors are represented by exactly the same numbers in the computer and we can safely argue, that to this accuracy the orbit is chaotic. The initial deviation vectors used for the orbit of Fig. 17c are and . From Eq. (17) it is evident that in the case of a Hamiltonian system with = 2 degrees of freedom, indices should decrease exponentially approaching the zero-value. In Fig. 17d we observe the behavior of for the same chaotic 2D orbit. Indeed, the indices obey to the exponential law given in Eq. (17). The evolution of the FNVI for this chaotic 2D orbit, for a time interval of time units, is presented in Fig. 3b, while dFNVI = 0.45. We see, that all the outcomes obtained using five different dynamical indicators coincide to the chaotic nature of the orbit.

We shall now proceed, presenting and comparing two more examples regarding 3D orbits of the Hamiltonian system (6). In Figure 18a-d, we can see the results provided by the four indicators for a regular 3D orbit with initial conditions: , , while the initial value of is found from the energy integral (6). The value of the energy is . In Fig. 18a we present the evolution of the DFLI for a time period of 5 time units. We observe that the DFLI fluctuates around small values (DFLI 10), which indicates regular motion. The initial deviation vector, , used for the orbit shown in Fig. 18a is . Fig. 18b depicts the evolution of the RLI for a time interval of time units, for the same regular orbit. The RLI clearly reveals the ordered character of the orbit, since RLI throughout the time evolution. Note that the axis is in log scale. The initial deviation vectors, used for the computation of the RLI are and , while . In Fig. 18c we present the evolution of SALI for a time interval of 5 time units. Also this method indicates the regular character of the orbit, since it exhibits small fluctuations around non-zero values (SALI ). The initial deviation vectors used for the orbit of Fig. 18c are and , but in general any other choice of initial values leads to similar behavior of the SALI. From Eqs. (18) and (19) it follows that in the case of a Hamiltonian system with = 3 degrees of freedom and will always remain different from zero, while , and should decay to zero following a power law, whose exponent depends on the number of deviation vectors that are initially tangent to the torus on which the orbit lies. Now, for the regular orbit of the 3D Hamiltonian (6) and for a random choice of initial deviation vectors, we expect the indices to behave as

(21) |

Fig. 18d shows the behavior of for the same regular 3D orbit. Indeed, the indices obey to the power laws respectively given in Eq. (21). In Fig. 18d with red color we plot the power law , the power law is colored in green, while the power law is plotted with blue color. Note that the axis is in log scale. With a more closer look in Fig. 18d, we observe that there is a superposition of and as they evolve almost identically. The evolution of the FNVI for this regular 3D orbit, for a time interval of time units, is presented in Fig. 4a, while dFNVI = 0.006. We see, that once more, all the outcomes derived using five different dynamical indicators coincide to the ordered nature of the orbit.

Finally, Figure 19a-d is similar to Fig. 18a-d but for a chaotic 3D orbit. In Fig. 19a-d, we present the results provided by the four dynamical indicators for a chaotic 3D orbit with initial conditions: , while the initial value of is found from the energy integral (6). The value of the energy is . In Fig. 19a we present the evolution of the DFLI for a time period of 5 time units. We observe that the DFLI displays a very rapid growth with time, which indicates chaotic motion. The initial deviation vector used for the orbit shown in Fig. 19a is . Fig. 19b depicts the evolution of the RLI for a time interval of time units, for the same chaotic orbit. Approximately for the first 200 time units, the RLI has low values (RLI ) but for the rest time interval RLI . Note that the axis is in log scale. Therefore, we may conclude that the orbit is chaotic. Note, that again the RLI method gives inconclusive or even misleading results for very small time intervals of the numerical integration. The initial deviation vectors, used for the computation of the RLI are and , while . In Fig. 19c we present the evolution of the SALI. We see that the SALI decreases very rapidly approaching the zero-value. At time units the SALI becomes zero, as it has reached the limit of the accuracy of the computer , which means that the two deviation vectors have the same direction. Thus, after time units the two normalized vectors are represented by exactly the same numbers in the computer and we can safely argue, that to this accuracy the orbit is chaotic. The initial deviation vectors used for the orbit of Fig. 19c are and . From Eq. (17) it is evident that in the case of a Hamiltonian system with = 3 degrees of freedom, all indices should decrease exponentially approaching the zero-value. In Fig. 19d we observe the behavior of for the same chaotic 3D orbit. Indeed, all the indices obey to the theoretically predicted law given in Eq. (17). The evolution of the FNVI for this chaotic 3D orbit, for a time interval of time units, is presented in Fig. 4b, while dFNVI = 0.52. We see, that all the outcomes obtained using five different dynamical indicators coincide to the chaotic nature of the orbit.

From the above examples, it becomes evident that all the different methods we used (DFLI, RLI, SALI and GALI) need not only the computation of the basic equations of motion but also the simultaneous computation of the so-called variational equations and of course several sets of deviation vectors. As the total number of the equations and inevitably the total number of the variables increases, we have also an increase to the time needed from the computational code in order to integrate the system of the differential equations and provide the results. Thus, our proposed FNVI method has one important advantage over all the other methods mentioned above, since it requires only the computation of the basic set of the equations of motion without any use of variational equations and deviation vectors. Let us present a specific illuminating example regarding the speed of the methods. Figure 20 presents the CPU time in h for each dynamical method required, using the same technique described in section 3 and of course the same processor, in order to construct the grid-plot shown in Fig. 7a. We observe, that the dFNVI method, needs only about 7 h of CPU time. The total time required by the spectrum is similar to the time needed for the dFNVI method. On the other hand, if we choose to use methods such as the DFLI, SALI and GALI, that need the additional computation of the variational equations, the total time required for the completion of the computational task increases significantly. Furthermore, we see from Fig. 20 that the most time consuming methods are the LCE and RLI. This is true, because when we use these two methods, we have to integrate each orbit at least for time units in order to obtain reliable information about the ordered or chaotic nature of the orbit. Thus, it becomes more than obvious, than the use of these two methods doubles the total CPU time needed by the dFNVI method. Therefore, our new method is not only reliable but also very fast compared with other methods and can be applied easily when we need to integrate a large number of orbits so as to construct a grid-plot of initial conditions (see Figs. 7 and 10).

Apart from the dynamical indicators described and applied previously, there are also other methods for distinguishing between ordered and chaotic motion, such as the mean exponential growth of nearby orbits (MEGNO) [4, 5], the power spectra of deviation vectors [30], the method of the low frequency power (LFP) [13, 29], the “0–1” test [10], as well as some other more recently introduced techniques [11, 24]. Moreover, a method that have also the ability to identify sticky orbits is the FtDt method [12]. This method uses the Fast Fourier Transform of a series of time interval, each one representing the time that elapsed between two successive points on the Poincaré surface of section. However, we do not feel that it is necessary to provide examples using all of these methods. We believe, that the comparison between the FNVI method and four well-known dynamical indicators proved its reliability and also revealed its remarkable efficiency.

## 6 Discussion and conclusions

In the present paper, we have tried to introduce a new, fast, efficient and easy to compute method in order to distinguish between order and chaos in 2D and 3D autonomous Hamiltonian systems. We have also conducted a detailed study of the behavior of this new indicator for both chaotic and regular orbits. The main results of this research can be summarized as follows:

1. The FNVI method proves to be an ideal detector of chaoticity independent of the dimensions of the dynamical system. It displays large, abrupt and random fluctuations for chaotic orbits, while it remains almost constant for ordered ones and so it clearly distinguishes between these two cases. Its main advantages are its simplicity, efficiency and reliability as it can rapidly and accurately determine the chaotic versus ordered nature of a given orbit. Of course, the FNVI method can provide only qualitative results regarding the nature of an orbit and therefore, we have to inspect the shape of FNVI each time in order to characterize an orbit. Obviously, this is not very practical when someone wants to check a large volume of orbits, so as to form an idea about the global structure of the dynamical system. Therefore, we proceeded one step further establishing a numerical criterion in order to quantify the results obtained by the FNVI method. This criterion derived by exploiting the shape of the FNVI. When the orbit is regular the FNVI remains almost constant, while in the case of a chaotic orbit it displays high fluctuations. Thus, we calculate the maximum and the minimum value of FNVI when . We choose this particular time interval because in the case where the orbit is regular, for time units there is a transient period of fluctuation, which may cause a problem to our criterion. Using the above procedure we defined the dFNVI. The value of the dFNVI can provide us the quantitative criterion that we seek. We point out, that it is not very easy to define a threshold value, so that the dFNVI being larger than this value reliably signifies chaoticity. Nevertheless, extensive numerical experiments in both dynamical systems indicate that in general, a good guess for this value could be 0.05. Thus, when dFNVI 0.05 the orbit is chaotic, while when dFNVI 0.05 the orbit ir ordered. This threshold value applies in both 2D and 3D dynamical systems.

2. We emphasize, that the main advantage of the FNVI method, is that it needs only the computation of the basic equations of motion. In particular, we only have to follow the evolution of the norm of the vector and compare its value in every time step throughout the entire time interval of the numerical integration with the initial value according to Eq. (4). The FNVI method is very fast, as it requires only time units of numerical integration, in order to provide reliable and conclusive results regarding the character of an orbit. Of course, there are also other methods and indicators capable to discriminate between ordered and chaotic orbits. However, the vast majority of them need the computation of the variational equations and inevitably the use of several sets of initial deviation vectors in order to function and provide their results. The additional use and therefore the computation of more equations increases significantly the total time needed from the computational routine to provide the output in each case (see Fig. 20). Thus, our proposed FNVI method has an important advantage over most the other dynamical methods, since it requires only the computation of the basic set of the equations of motion. Therefore, our method is very fast which makes it is an ideal tool when we need to calculate a large number of orbits so as to construct a grid of initial conditions.

3. Using the numerical quantified version of FNVI, that is the dFNVI, we are in a position to characterize reliable an orbit of being chaotic or ordered. Exploiting the advantages of the dFNVI method, we have constructed detailed phase-space portraits (grid-plots) both for 2D and 3D Hamiltonian systems, where the chaotic and ordered regions are clearly distinguished (see Figs. 7 and 10). We were thus able to trace in a fast and systematic way very small islands of ordered motion, whose detection by traditional methods would be very difficult, if not impossible, and moreover very time consuming. This approach is therefore, expected to provide useful tools for the location of stable periodic orbits, or the computation of the phase-space volume occupied by ordered or chaotic motion in multi-dimensional systems ), where the PSS plots can not be easily visualized and interpreted and furthermore, very few other similar techniques of practical value are available. Here, we have to point out that grid-plots can be also constructed using other dynamical indicators. For the reasons mentioned in the previous point, the dFNVI method is much more faster than most of the other methods.

4. We used the new FNVI method in order to follow the time evolution of a sticky orbit. The study of the evolution of sticky orbits (orbits that, although they are chaotic, they are restricted in a thin layer of phase space for a large value of integration time) is crucial for the understanding of the structure of the dynamical system. The ability of our new method to extract reliable information concerning the chaotic or regular behavior of an orbit enables us to overcome the ambiguity in the case of sticky orbits. In the case of sticky orbits, the distinction can not be easily made using the PSS technique, as a sticky orbit can be mistaken for an ordered one, if not enough time is given and the orbit has not yet escaped to the surrounding chaotic domain (see Fig. 13a-d). Moreover, the LCE is not sensitive enough to identify the difference between a regular and a sticky orbit (see Fig. 15a). On the contrary, the FNVI method can be applied to obtain a clear and reliable distinction between ordered and sticky orbits (see Fig. 15b) and this could be regarded as a major advantage of our method.

In the present research, we have tested approximately orbits (2D and 3D) in the dynamical systems (6) and (7). Since our results are consistent, we may safely conclude that the FNVI method (and also its numerical version the dFNVI) is a very fast and reliable tool for distinguishing between order and chaos in Hamiltonian systems. It is our future plans, to apply this new method in other kinds of interesting potentials and also in more complicated dynamical systems, in order to test further its efficiency and reliability. Furthermore, we will apply this new method in order to study the evolution of three-dimensional sticky orbits. Moreover, with additional theoretical work, we will investigate if it always discriminates correctly between ordered and chaotic motion in Hamiltonian systems, especially in dynamical systems with three degrees of freedom.

## Acknowledgements

I would like to express my warmest thanks to Dr. D.D. Carpintero for all the illuminating and creative discussions during this research. The author would also like to thank the two anonymous referees for the careful reading of the manuscript and for their very useful and aptly suggestions and comments which allowed us to improve significantly the quality and the clarity of the present article.

### Footnotes

- journal: Nonlinear Dynamics

### References

- Arribas, M., Elipe, A., Floria, L., Riaguas, A.: Oscillators in resonance . Chaos, Solitons & Fractals 27, 1220-1228 (2006)
- Benettin, G., Galgani, L., Giorgilli, A., Strelcyn, J.-M.: Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems - A method for computing all of them. I - Theory. II - Numerical application. Meccanica 15, 9-30 (1980)
- Caranicolas, N.D., Zotos, E.E.: Investigating the nature of motion in 3D perturbed elliptic oscillators displaying exact periodic orbits. Nonlinear Dyn. 69, 17951805 (2012)
- Cincotta, P.M., Simó, C.: Simple tools to study global dynamics in non-axisymmetric galactic potentials - I. A&A Supplement Series 147, 205-228 (2000)
- Cincotta, P.M., Giordano, C.M., Simó, C.: Phase space structure of multi-dimensional systems by means of the mean exponential growth factor of nearby orbits. Physica D 182, 151-178 (2003)
- Contopoulos, G., Grousousakou, E., Voglis, N.: Invariant spectra in Hamiltonian systems. A&A 304, 374-380 (1995)
- Elipe, A., Deprit, A.: Oscillators in resonance. Mech. Res. Commun. 26, 635-640 (1999)
- Froeschlé, C., Lega, E., Gonczi, R.: Fast Lyapunov Indicators. Application to Asteroidal Motion. Celest. Mech. Dyn. Astron. 67, 41-62 (1997)
- Froeschlé, C., Gonczi, R., Lega, E.: The fast Lyapunov indicator: a simple tool to detect weak chaos. Application to the structure of the main asteroidal belt. P&SS 45, 881-886 (1997)
- Gottwald, G.A., Melbourne, I.: A new test for chaos in deterministic systems. RSPSA 460, 603-611 (2004)
- Howard, J.E.: Discrete Virial Theorem. Celest. Mech. Dynam. Astron. 92, 219-241 (2005)
- Karanis, G.I., Vozikis, Ch.L.: Fast detection of chaotic behavior in galactic potentials. Astron. Nachr. 329, 403-412 (2008)
- Kotoulas, T., Voyatzis, G.: Comparative Study of the 2:3 and 3:4 Resonant Motion with Neptune: An Application of Symplectic Mappings and Low Frequency Analysis. Celest. Mech. Dynam. Astron. 88, 343-363 (2004)
- Laskar, J.: The chaotic motion of the solar system - A numerical estimate of the size of the chaotic zones. 1990 Icarus 88, 266-291 (1990)
- Laskar, J.: Hamiltonian Systems with Three or More Degrees of Freedom. Plenum, New York p 134 (1999)
- Laskar, J., Froeschlé, C., Celletti, A.: The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping. Physica D 56, 253-269 (1992)
- Lieberman, M.A., Lichtenberg, A.J.: Regular and Chaotic Dynamics. Springer, Berlin (1992)
- Meiss, J.D.: Symplectic maps, variational principles, and transport. RvMP 64, 795-848 (1992)
- Papaphilippou, Y., Laskar, J.: Frequency map analysis and global dynamics in a galactic potential with two degrees of freedom. A&A 307, 427-449 (1996)
- Papaphilippou, Y., Laskar, J.: Global dynamics of triaxial galactic models through frequency map analysis. A&A 329, 451-481 (1998)
- Patsis, P.A., Efthymiopoulos, C., Contopoulos, G., Voglis, N.: Dynamical spectra of barred galaxies. A&A 326, 493-500 (1997)
- Saito, N., Ichimura, A.: In: Casati, G., Ford, J. (eds.) Stochastic Behavior in Classical and Quantum Hamiltonian Systems. Lecture Notes in Physics, vol. 93, p. 137. Springer, Berlin (1979)
- Sándor, Zs., Érdi, B., Széll, A., Funk, B.: The Relative Lyapunov Indicator: An Efficient Method of Chaos Detection. Celest. Mech. Dyn. Astron. 90, 127-138 (2004)
- Sideris, I.V.: In Gottesman, S.T., Buchler, J-R., et al. (Eds.), Nonlinear Dynamics in Astronomy and Astrophysics, In: Annals of the New York Academy of Science vol. 1045 The New York Academy of Sciences, New York (2005)
- Skokos, Ch.: Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits. J. Phys. A 34, 1002910043 (2001)
- Skokos, Ch., Antonopoulos, Ch., Bountis, T.C., Vrahatis, M.N.: Detecting order and chaos in Hamiltonian systems by the SALI method. J. Phys. A 37, 6269-6284 (2004)
- Skokos, Ch., Bountis, T.C., Antonopoulos, Ch.: Geometrical properties of local dynamics in Hamiltonian systems: the generalized alignment index (GALI) method. Physica D 231, 30-54 (2007)
- Voglis, N., Contopoulos, G.: Invariant spectra of orbits in dynamical systems. J. Phys. A 27, 4899-4912 (1994)
- Voyatzis, G., Ichtiaroglou, S.: On the spectral analysis of trajectories in near-integrable Hamiltonian systems. J. Phys. A 25, 5931-5943 (1992)
- Vozikis, Ch.L, Varvoglis, H., Tsiganis, K.: The power spectrum of geodesic divergences as an early detector of chaotic motion. A&A 359, 386-396 (2000)
- Vrahatis, M.N., Bountis, T.C., Kollmann, M.: Periodic Orbits and Invariant Surfaces in 4-D Nonlinear Mappings. Int. J. Bifurc. Chaos 6 1425-1437 (1996)
- Vrahatis, M.N., Isliker, H., Bountis T.C.: Structure and Breakdown of Invariant Tori in a 4-D Mapping Model of Accelerator Dynamics. Int. J. Bifurc. Chaos 7 2707-2722 (1997)
- Zotos, E.E.: A new dynamical model for the study of galactic structure. New Astronomy 16, 391-401 (2011)
- Zotos, E.E.: A New Dynamical Parameter for the Study of Sticky Orbits in a 3D Galactic Model Baltic Astronomy 20, 339-354 (2011)
- Zotos, E.E.: Application of new dynamical spectra of orbits in Hamiltonian systems. Nonlinear Dyn. 69, 20412063 (2012)