Stickiness and recurrence plots: An entropy-based approach

The phase space of a typical quasi-integrable Hamiltonian system is, in general, neither integrable nor uniformly hyperbolic, but there is a coexistence of chaotic and regular domains.³ The regular dynamics consists of periodic and quasiperiodic orbits that lie on invariant tori, while the chaotic orbits fill densely the available domain in phase space. For the special case of two-dimensional area-preserving maps, the existence of islands of regularity filled with Kolmogorov–Arnold–Moser (KAM) invariant tori separates the phase space into distinct regions; i.e., orbits in the chaotic sea will never enter any island, and the periodic and quasiperiodic orbits inside of an island will never reach the chaotic sea.^3,4 Due to the existence of islands embedded in the chaotic sea, the latter constitutes a fat fractal⁵ and it is challenging to determine exactly the islands’ boundary. The islands are surrounded by smaller islands, which are in turn surrounded by even smaller islands. This structure repeats itself for arbitrarily small scales giving rise to an infinite hierarchical islands-around-islands structure.⁶ 

The phenomenon of stickiness^7–11 emerges due to this complex interplay between islands and chaotic sea. The chaotic orbits that approach an island may spend an arbitrarily long, but finite, time in its neighborhood, in which the orbits will behave similarly as quasiperiodic orbits. Before escaping, the orbits are trapped in a region bounded by cantori,^4,10,12 which are a Cantor set, formed by the remnants of the destroyed KAM tori. Unlike the KAM tori that are full barriers to the transport in phase space, the cantori act as partial barriers, where the orbits may be trapped in the region bounded by them, and once trapped inside a cantorus, the chaotic orbits may cross an inner cantorus and so on to arbitrarily small levels in the hierarchical structure of islands-around-islands.

Since in two-dimensional area-preserving maps, the islands and the chaotic sea are distinct and disconnected domains, it is of major importance to characterize orbits into these two categories. The traditional and most known method to characterize the dynamics of a system is through the evaluation of the Lyapunov exponents.^13,14 For a two-dimensional mapping, there are two Lyapunov exponents, and ⁠, and the dynamics is chaotic if one of them is positive. When the mapping is area-preserving, the sum of the two exponents must be zero; i.e., ⁠. In this scenario, the regular orbits have Lyapunov exponents equal to zero for infinite times, while the chaotic ones exhibit ⁠.

If sticky regions are present in the phase space, the Lyapunov exponents may not be the optimal choice to detect chaotic orbits due to the trappings around the islands. When the orbit is trapped, the largest Lyapunov exponent decreases and this makes its convergence slower; i.e., it takes longer to reach the asymptotic (infinite-time) value. As an alternative, a new method based on ergodic theory has recently been proposed to detect chaotic orbits.^15–18 It proved to be a better option to distinguish between regular and chaotic orbits than the Lyapunov exponents. However, both this new method and the Lyapunov exponents require very long time series in order to obtain reliable accuracy. When only a short time series is available, a possible approach is to use the recurrence quantification analysis (RQA).^19–23 The RQA was developed to quantify the dynamics of a system by means of the recurrences of the orbit in phase space.

The most used RQA measures (e.g., the recurrence rate and the determinism) can, in some sense, detect different transitions occurring in nonlinear systems. However, we seek a measure based on an intrinsic property of dynamical systems: quasiperiodic orbits lying on invariant circles can have at most three different return (recurrence) times, which is the time needed for the orbit return to a given neighborhood of a point in the orbit, as stated by Slater’s theorem.^1,2 We can obtain the recurrence times by simply defining a recurrence region and counting how long it takes to the orbit to return to this given region. It is also possible to use the recurrence plots (RPs) to estimate the recurrence times: the white vertical lines in the RP give us a lower estimate of the recurrence times.^24–28 Thus, in this paper, we propose using RPs to characterize the dynamics of a nonlinear system. We focus on a further type of RP-based measure to identify regular and chaotic regions and sticky regions as well, namely, the Shannon entropy of the recurrence times: the recurrence time entropy (RTE).^29,30 The RTE was originally introduced without any connection to the RPs,²⁹ and it was shown that it can provide a good estimate for the Kolmogorov–Sinai entropy²⁷ and the largest Lyapunov exponent.³¹ 

We find that, with the RTE, we can identify very clearly the regular regions and the transitions to chaotic motion as one parameter of the system is varied. We also find that, by computing the finite-time RTE distribution, we can identify a multi-modal distribution, in which each maximum is related to a different hierarchical level in the islands-around-islands structure. Moreover, each of these regions corresponds to a power-law decay of the cumulative distribution of trapping times.

This paper is organized as follows. In Sec. II, we introduce the standard map and briefly comment on some of the properties of two-dimensional quasi-integrable Hamiltonian systems. We also discuss a few approaches of how to detect sticky orbits in the phase space of such systems. In Sec. III, we introduce the concept of recurrence plots and show that it can be used to characterize the dynamics of a dynamical system. In this section, motivated by Slater’s theorem, we also propose the use of the RTE, estimated from the RP, to quantify the dynamics as periodic, quasiperiodic, and chaotic. In Sec. IV, we apply this entropy-based measure for the standard map and show that it is positively correlated to the largest Lyapunov exponent, with a high correlation coefficient, and also show that it is possible to detect and characterize different regimes of stickiness present in the dynamics. Section V contains our final remarks.

In spite of its simple mathematical form, the standard map exhibits all the features of a typical quasi-integrable Hamiltonian system, and it has become a paradigmatic model for the study of properties of chaotic motion in quasi-integrable Hamiltonian systems. For ⁠, the dynamics is regular, the system is integrable, and all orbits lie on invariant rotational tori. As the nonlinearity parameter increases, the “sufficiently irrational” invariant tori persist, as predicted by the KAM theorem,³ whereas the rational tori are destroyed. When is sufficiently large, it turns out that all the invariant rotational tori are destroyed. For the standard map, the last invariant rotational torus ceases to exist for the critical value ⁠,³³ leading to a scenario of global stochasticity.

One of the main features of quasi-integrable Hamiltonian systems is the phenomenon of stickiness. The phase space of a quasiperiodic orbit (blue), a chaotic orbit (black), and a sticky orbit (red) of the standard map (1) with and the initial conditions ⁠, and ⁠, respectively, are displayed in Fig. 1. The orbits were iterated for times, and we also calculated the largest Lyapunov exponent, ⁠, of each orbit, namely, ⁠, ⁠, and ⁠, respectively. For this value of ⁠, one single chaotic orbit fills a significant portion of phase space (black dots), and the most prominent sticky region is around the period-6 satellite islands (red dots). An orbit initialized in this region is trapped for a long time until it escapes to the chaotic sea, and this affects some properties of the orbit, the ⁠, in particular. Although chaotic and sticky orbits both have ⁠, the sticky orbit has a lower value of than the chaotic one. The for the quasiperiodic orbit is small, but not exactly zero due to the finite iteration time ⁠: as ⁠, ⁠. Furthermore, the closer the quasiperiodic orbit is to the elliptic point, the faster the convergence of toward zero.³⁴ 

The stickiness is usually characterized through the recurrence time statistics,^6,11,35–40 although other methods have been proposed, such as the finite-time Lyapunov exponent (FTLE),^41–44 the rotation number,⁴⁵ and the recurrence quantification analysis (RQA).^24–26,46 Szezech et al.⁴¹ showed that for the case where the phase space has stickiness regions, the distribution of the finite-time Lyapunov exponents is bimodal. The finite-time Lyapunov exponent has also been used to characterize stickiness in high-dimensional Hamiltonian systems.^43,44 More recently, Santos et al.⁴⁵ showed that the rotation number is a faster method, compared to the finite-time Lyapunov exponent, to verify the presence of sticky orbits in the phase space. The most commonly used measures of the RQA, such as the determinism and the recurrence rate, can also characterize stickiness in a similar way.^24,25,46

In Sec. III, we review the concept of recurrence plots, and we propose a quantity based on them, different from the traditional measures of RQA, to characterize the dynamics of a nonlinear system.

The recurrence matrix is a symmetric, binary matrix that contains the value for recurrent states and the value for nonrecurrent ones. Two states are said to be recurrent when the state at is “close” (up to a distance ⁠) to a different state at ⁠, i.e., ⁠. The choice of the threshold is not arbitrary. If is chosen too large, almost every point is recurrent with every other point. On the other hand, if is chosen too small, there will be almost no recurrent states. Several rules have been proposed. Some consider with a fixed recurrence point density of the RP.⁴⁸ Another possibility is to consider as a fraction of the standard deviation, ⁠, of the time series.^49,50 Regardless of the choice we make, the effect of a finite will never disappear; a new study has shown that using is not the best choice,⁵¹ and in this work, we consider the threshold to be of the time series standard deviation, considering the maximum norm approach (see the  Appendix); i.e., ⁠. For a detailed discussion about the effect of in our results, see the  Appendix.

Graphically, the recurrent states are represented by a colored dot, and the recurrence matrix displays different patterns according to the dynamics of the underlying system. In Fig. 2, we show the recurrence matrix for the first iterations of (a) a quasiperiodic orbit, (b) a chaotic orbit, and (c) a sticky orbit of the standard map with shown in Fig. 1. The RP of the quasiperiodic orbit consists mainly of uninterrupted diagonal lines [Fig. 2(a)]. The vertical distance between these lines is regular and corresponds to the different return times of the orbit,^24–26,28 whereas the RP of the chaotic orbit displays short diagonal lines and the vertical distances between them are not as regular as in the quasiperiodic case. The RP of the sticky orbit seems to be between these two cases. The diagonal lines are longer than those in the chaotic case, indicating that the sticky orbit is more regular than the chaotic one, but the diagonal lines are not as long as those in the quasiperiodic case. Moreover, the vertical distances between the diagonal lines have some regularity.

Therefore, we can use the RP to quantify the dynamics of the system. Several measures based on the length of the diagonal lines and based on the length of the vertical lines have been proposed, such as the recurrence rate, the determinism, the laminarity, the maximal length of diagonal and vertical lines, among others. A complete discussion about these and other quantifiers can be found in Refs. 19–23 and references therein.

In this way, a periodic orbit, which has only one return time (the period itself), will have ⁠. A quasiperiodic orbit, which has three return times, will lead to a low value of ⁠, whereas a chaotic orbit will be characterized by a high value of ⁠. As has been stated, the chaotic orbit that experiences stickiness spends an arbitrarily long time in the neighborhood of an island in which the orbit exhibits similar behavior as a quasiperiodic orbit. This fact can be seen from the recurrence matrix shown in Fig. 2(c). Hence, we expect the RTE of sticky orbits to be smaller than the chaotic ones, but higher than it would be for a quasiperiodic orbit.

In this section, we evaluate the for the standard map, and we show that the RTE can be used to characterize the dynamics of the system and to detect the presence of stickiness regions. Unless mentioned otherwise, we use a time series of size for all our simulations.

In Fig. 3, we show the and the RTE of the standard map for a fixed initial condition as a function of the nonlinearity parameter ⁠. We notice that whenever ⁠, the RTE is large. Also, we can identify the windows of regularity, in which is zero and the RTE assumes low values. Furthermore, even though the Lyapunov exponent goes faster to zero around the elliptic point,³⁴ only with the values of the Lyapunov exponent, we cannot distinguish between the periodic and quasiperiodic orbits in two-dimensional Hamiltonian systems. There are several values of for which ⁠, indicating that at these points, the initial condition is very close to a periodic orbit [Fig. 3(b)].

To have a better visualization of the correspondence between the and the ⁠, we plot in Fig. 4 the values of and for a grid of initial conditions uniformly distributed in the phase space with and in the parameter space with ⁠. Figures 4(b) and 4(e) are a magnification of one of the period-6 satellite islands of Figs. 4(a) and 4(d), respectively. We see that the RTE captures all the features of the Lyapunov exponent but even more. In the chaotic sea, where is large, is also large and inside the islands, where ⁠, the RTE is low. In addition to that, in the regions where the rotation number of an orbit is close to a rational number, the RTE is smaller (blue to purple) [Fig. 4(d)].^15,61 The RTE decreases as we get closer to the elliptic point, as we can see in Figs. 4(e) and 4(f). In the latter, we can see the transitions from regular to chaotic behavior, where bifurcations occur as changes.

For the chosen parameter values of the standard map, ⁠, it is known that the system exhibits the stickiness effect, and due to that, the distribution of the finite-time Lyapunov exponent is bimodal.⁴¹ The most prominent sticky region for this parameter is the region in between the main island and the period-6 satellite islands. The decreases in this region when compared with the rest of the chaotic sea [Figs. 4(a) and 4(b)]. With the RTE, we observe the same behavior [Figs. 4(d) and 4(e)].

Next, we consider a single chaotic orbit of the standard map, and we follow up its evolution for a long iteration time ⁠. For long times, the trajectory fills the entire chaotic component of the phase space, and usually, the stickiness acts for very long, but finite, times before an orbit escapes to the chaotic sea. Thus, the transitions from different regimes in the dynamics of the orbit can be better understood considering “finite-time” ⁠. Therefore, we compute the RTE along the evolution of a single chaotic orbit in windows of size ⁠, ⁠, where ⁠, and define the probability distribution of the finite-time RTE, ⁠, by computing a frequency histogram of such as the one shown in Fig. 5(a) for and ⁠. To compute this distribution, we use the matplotlib function hist^62,63 with the parameter bins=”auto”. It uses the maximum of Sturge’s rule and the Freedman–Diaconis estimator to compute the number of bins taking into account data variability and data size.⁶⁴ The inset in Fig. 5(a) shows the finite-time RTE “time series” for the interval from to ⁠. We see abrupt changes in the value of ⁠, indicating the transitions from different regimes in the dynamics of the orbit. These changes in the value of the RTE cause its probability distribution to split into more than one mode. Szezech et al.⁴¹ reported that for this value of ⁠, the distribution of the finite-time Lyapunov exponent is bimodal. What we see is that the minor peak of Fig. 3 in Ref. 41 consists, in fact, of multiple peaks, as suggested by Harle and Feudel.⁴² 

The multi-modal distribution is due to the infinite hierarchical islands-around-islands structure embedded in the phase space. When the orbit is in the chaotic sea, the time- RTE is high, corresponding to the largest maximum of the distribution. On the other hand, when the orbit is trapped near an island, the RTE is low and the distribution exhibits smaller maxima for small values of ⁠. Once trapped in the neighborhood of an island, the orbit may enter an inner level in the hierarchical structure, and these transitions to different levels are the cause of the multi-modal distribution.⁴² Moreover, the closer to zero is ⁠, the higher the hierarchical level of the island on which neighborhood the orbit is trapped.⁶⁵ Hence, multiple peaks are formed for small values of RTE [Fig. 5(a)].

In order to identify the regions in phase space that correspond to the peaks in the distribution, we monitor the time series and plot the phase space positions with different colors for different ranges of ⁠. The blue, red, green, and black points represent the phase space position when ⁠, ⁠, ⁠, and ⁠, respectively [Figs. 5(b) and 5(c)]. Different peaks are indeed related to different hierarchical levels of the structure, and by means of the RTE, it is possible to distinguish between them very well. Additionally, the black and green points shadow the manifolds along which the non-trapped orbits leave the sticky region.⁶⁶ Even though only the corresponding hierarchical levels around the period-6 island are shown [Fig. 5(c)], all island chains contribute to the finite-time RTE distribution. For completeness, we also investigate the phase space points that generate the peak for high values of RTE. Whenever ⁠, we plot the phase space points [Fig. 6(a)]. The phase space components in Figs. 5(b) and 6(a) complement each other: in Fig. 6(a), the chaotic sea is shown, i.e., the hyperbolic component of phase space, whereas in Fig. 5(b), we see the nonhyperbolic component. Nonhyperbolicity can inhibit chaotic orbits from visiting some regions due to tangencies between stable and unstable manifolds.^67–69 

Several approaches to detect the existence of sticky orbits in the phase space of two-dimensional area-preserving systems have already been proposed and studied in previous studies. RP-based measures have also been used for this purpose, however, without as much care as the Lyapunov exponents, for example. In this paper, we have proposed the use of the Shannon entropy of the distribution of the recurrence times, RTE, estimated from the recurrence plots, to detect and characterize the stickiness effect in the standard map. We have shown that the RTE is positively correlated to the largest Lyapunov exponent, with a high correlation coefficient, and that it is possible to distinguish among the different types of motion using this entropy-based measure.

It is well-known that for systems that exhibit the stickiness effect, the transitions from fully chaotic motion to different levels in the hierarchical structure of islands-around-islands cause the FTLE distribution to have more than one peak. In fact, the distribution is multi-modal, in which each peak represents a different hierarchical level of islands. For the chosen nonlinearity parameter, ⁠, it was previously reported that the FTLE is bimodal.⁴¹ However, we have shown here that the maximum for small values of / is, in fact, composed of several minor peaks, as suggested by Harle and Feudel.⁴² This suggests that the RTE can be an alternative way of characterizing the stickiness effect, with which it is possible to distinguish among the different hierarchical levels in the islands-around-islands structure embedded in the chaotic component of phase space.

The presence of sticky regions in phase space affects global properties of the system, such as the distribution of the trapping times. We have shown that, by monitoring the finite-time RTE time series, it is possible to collect a set of trapping times for each of the different hierarchical levels detected by the finite-time RTE distribution and obtain the probability distribution of ⁠, ⁠, and its cumulative distribution, ⁠, of each hierarchical level. of fully chaotic systems has an exponential decay, whereas exhibits a power-law tail when sticky regions are present in the phase space. After separating among the distinct hierarchical levels, we have shown that the cumulative distribution of the hyperbolic component has indeed an exponential decay and that the power-law tail, characteristic of stickiness, is observed in the distribution of the different hierarchical levels.

One interesting point we plan to investigate in the future is whether the RTE can characterize also higher dimensional systems (e.g., 4D symplectic maps) using the methodology presented in this paper. Another interesting study is the critical value ⁠,³³ where the last invariant rotational torus ceases to exist. For this parameter, there is a sequence of cantori, and this could give an interesting histogram for the finite-time RTE.

All of our simulations regarding the evaluation of the recurrence matrix were made using the pyunicorn package,⁷⁰ and even though these results were obtained with the standard map, we expect similar results for any quasi-integrable Hamiltonian system that exhibits stickiness.

We wish to acknowledge the support of the Araucária Foundation, the Coordination of Superior Level Staff Improvement (CAPES), under Grant No. 88881.143103/2017-1; the National Council for Scientific and Technological Development (CNPq), under Grant Nos. 403120/2021-7, 301019/2019-3; and the São Paulo Research Foundation (FAPESP) under Grant Nos. 2018/03211–6 and 2022/04251-7. J.K. has been supported by the Alexander von Humboldt Polish Honorary Research Scholarship 2020 of the Foundation for Polish Science. We would also like to thank the 105 Group Science⁷¹ for fruitful discussions.

The authors have no conflicts to disclose.

Matheus R. Sales: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Michele Mugnaine: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). José D. Szezech Jr.: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Ricardo L. Viana: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Iberê L. Caldas: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Norbert Marwan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jürgen Kurths: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The source code to reproduce the results of this paper is freely available in the Zenodo archive.⁷² 

In Sec. III, we introduced the concept of RPs and made some considerations regarding the choice of the threshold ⁠: we chose to be of the time series standard deviation, ⁠. In the following, we provide an analysis of the effect of in our results.

To determine the optimal method for calculating and an appropriate value of ⁠, we compute the correlation coefficient, Eq. (6), between and RTE as a function of (in units of ⁠) using the concatenation and norm approaches to calculate [Fig. 7(b)]. The three approaches yield similar correlation coefficients. Even if we choose a very small (⁠ of ⁠), we still obtain a high correlation coefficient (⁠⁠). Also, increasing does not affect significantly ⁠, indicating that in our case, the choice of is not that sensible as it would be in other cases. In fact, there is a range of values for in which the results are good.

Even though the three approaches of calculating the standard deviation give similar results in our case, choosing the concatenation approach when one time series has a different value range than the others might strongly bias the standard deviation. Therefore, in our opinion, the most appropriate approach is the norm approach (either Euclidean or maximum).