Hidden Strange Nonchaotic Attractors

In this paper, it is found numerically that the previously found hidden chaotic attractors of the Rabinovich–Fabrikant system actually present the characteristics of strange nonchaotic attractors. For a range of the bifurcation parameter, the hidden attractor is manifestly fractal with aperiodic dynamics, and even the finite-time largest Lyapunov exponent, a measure of trajectory separation with nearby initial conditions, is negative. To verify these characteristics numerically, the finite-time Lyapunov exponents, ‘0-1’ test, power spectra density, and recurrence plot are used. Beside the considered hidden strange nonchaotic attractor, a self-excited chaotic attractor and a quasiperiodic attractor of the Rabinovich–Fabrikant system are comparatively analyzed.


Introduction
In the last three decades, strange nonchaotic attractors (SNAs) have attracted much attention from both theoretical and experimental points of view, reportedly arising in many physically relevant situations. Though somewhat exotic, SNAs are not at all rare. In fact, SNAs have been reported in quasiperiodically forced pendulum, quantum particles, quasiperiodic potentials, biological oscillators, quasiperiodically driven Duffing-type oscillators, velocity dependent oscillators, electronic circuits, etc. (see, e.g., [1] and references therein). Additionally, some models are studied in higher dimensions and for both discrete and continuous systems. Commonly existing in quasiperiodically driven nonlinear systems, SNAs are found and studied for non autonomous quasiperiodically forced dynamical systems and can be regarded as "intermediate" structures linking regularity and chaos. However, as this paper shows, there are undriven continuous systems with SNAs, or even with hidden SNAs.
SNAs were introduced in 1984, when Grebogi, Ott, Pelikan, and Yorke constructed dynamical systems with attractors that were strange but not chaotic [2] (see also [3][4][5][6][7]). Results on flows containing SNA were obtained far before than the term SNA was coined (see [8], where some references (in Russian) can be found)). The first experimental observation of an SNA was in a magnetoelastic ribbon [9] (see also [10][11][12]). The 'strange' characteristic refers to the fractal, nontrivial, and complicated geometry of the attractor. There are authors who consider SNAs to be "strange" because they are not piecewise differentiable. The relatively new term, SNAs, has still not been precisely formulated and their existence has not been rigorously proven. Therefore, most results on SNAs, such as those in this paper, are based on numerical analysis.
The 'nonchaotic' property means that SNAs do not show sensitivity to initial conditions-i.e., with negative Lyapunov exponents (LEs), just like regular systems. Note that some authors consider the nonchaotic property to be defined by nonpositive LEs (see, e.g., [6]). Although it is difficult to distinguish an SNA visually from a strange chaotic attractor, dynamically there are important differences. Indeed, orbits on SNAs and chaotic attractors are non-periodic, but compared to chaotic attractors the trajectories on SNAs are not separate from each other. The dynamical behavior characterizing SNAs occurs in a very narrow range of the control parameter values.
In principle, strange nonchaotic attractors occur in all dissipative dynamical systems that exhibit the period-doubling route to chaos, where the attractors formed at the accumulation points are fractal sets with zero Lyapunov exponents. Such attractors are, however, not physically observable because the set of parameter values for them to arise has a Lebesgue measure of zero in the parameter space [13].
SNAs are typically not robust in the absence of quasiperiodic driving. Therefore, they are not expected to occur naturally, whilst, as shown by Grebogi et al. in [2], if the system is quasiperiodically driven then SNAs can be robust. A robust SNA is typical but not vice versa.
From a computational point of view, attractors can be classified as self-excited (SE) and hidden. While SE attractors can be localized numerically with some standard computational procedures by starting from a point in a neighborhood of unstable equilibrium, the basins of hidden attractors, which can be chaotic (HCAs) or non-chaotic (e.g., stable cycles), are not connected with equilibria. Thus, if the trajectories are attracted by some attractor, the attractor is an SE attractor, while a hidden attractor cannot be excited by the unstable equilibria (if they exist) and trajectories from neighborhoods of unstable equilibria are attracted either by some attractors (stable equilibria or cycles) or diverge to infinity [14][15][16][17][18][19]. Therefore, the HCAs are usually found empirically. In other words, it means to verify if the attraction basins of the HCA do not intersect small neighborhoods of unstable equilibria. However, finding the attraction basins of HCAs analytically remains an open and challenging problem.
The investigation of hidden oscillations actually arose in the second part of Hilbert's 16th problem and in Bautin's works, related to nested limit cycles in quadratic systems, which showed the necessity of studying hidden oscillations to solve this problem [15]. HCAs may appear in, e.g., chaotic systems with one or several stable equilibriums or with a line equilibria (infinitely many equilibria), or without equilibrium [20][21][22][23].
In this paper, attractors and transients of the Rabinovich-Fabrikant (RF) are differentiated via their lifetime and it is shown numerically that some HCAs of the RF system present the characteristics of SNAs. For this purpose, the spectrum Λ of the finite-time local Lyapunov exponents (LEs), the '0-1' test (Appendix A), the two-sided Power Spectral Density (PSD) determined for one component, the Correlation Dimension (D 2 ), and the Recurrence Plot (RP) (Appendix B) are used. Additionally, histograms and the difference between two close trajectories are utilized to differentiate SNAs from chaotic attractors. This paper is organized as follows: Section 2 presents the Rabinovich-Fabrikant system with some representative attractors and Section 3 shows that the system has hidden strange nonchaotic attractors (HSNAs). The Appendixes A-C present the notions of the '0-1' test and the recurrence plot. Finally, the conclusion in Section 4 ends the paper.

Hidden Attractors in RF System
The RF system [24,25], arising from the modulation instability in a nonequilibrium dissipative medium, is modeled by: where a and b are some real positive parameters and b is the bifurcation parameter.
The system has extremely rich dynamics, presenting coexisting attractors, SE chaotic attractors (SECAs), and HCAs [26]. The equilibria of the system are: and X * 0 (0, 0, 0), with: Due to the general instability of chaotic solutions and due to the particular sensibility of the RF system to numerical integration, the dedicated numerical methods implemented in known software cannot be used here with 100% confidence, especially for long-time integration, because global errors increase with the size of the integration interval [27,28]. To numerically integrate the RF system, in this paper Matlab solver ode45 is used. To note that the system presents some stiffness characteristics (the x 3 component of the solution decreases faster than x 1,2 [29]). However, the matlab solvers for stiff systems give, in this case, similar results). Most importantly, the numerical results cannot be improved by decreasing the integration step, since the integration error has an extremum as a function of the integration step-size [30]. Another approach would be higher-accuracy calculations, but this approach has restrictions: on one hand, the way to decrease the error is narrow; on another hand, the number of operations needed, if a very small integration step-size is used, could increase dramatically (see [30] and the references therein). Additionally, note that it is possible to have reliable numerical simulations only in a relatively small time interval [31][32][33]. For example, for the classical Lorenz system, it is a real challenge to get solutions with small errors (in the sense that the largest error between two computed solutions for small but different time-steps has a tolerance of 5 × 10 −14 ) for t ∈ [0, 100] [31].
Therefore, one can never be confident if an apparently chaotic numerical trajectory on larger time intervals represents a real chaotic trajectory or just a chaotic transient. Beside these obstacles in the numerical integration of ODEs of the RF system, the first intensive numerical investigations [25] revealed an extremely strong dependence of the solution on the utilized numerical methods, initial conditions, and the parameter b.
On the other hand, transient chaos appears when trajectories starting from some random initial conditions are chosen near the boundary of the basin of attraction (boundary crisis) or are very close to some bifurcation point. The underlying trajectories behave chaotically, possibly for a relative long time, and then quite abruptly switch to some final periodic state [34][35][36]. Because the lifetime of the transients can be very long (thousands and tens of thousands of numerical iterations), one can never be 100% confident that an observed chaotic system is not a transient (Sprott's communication). Therefore, to differentiate transients from attractors of the RF system numerically, a compromise with a sufficiently large time interval is made, with [0, T] and T = 10, 000. For a = 0.1, b ∈ [0.2808, 0.2907] and with the initial condition x 0 = [0.02, 0.12, 0.53], Figure 1 shows the lifetime of the trajectories together with the variation in the spectrum of the LEs, In order to increase T over T = 10, 000, options = odeset('RelTol , 1e − 6, 'AbsTol , 1e − 9) is used [37]. Note that without using options, some of attractors of the RF system transform into chaotic transients. Thus, in this paper the trajectories which reach a final regular motion before T = 10, 000 (here the stable equilibria X * 1,2 ) with an error of less than 1 × 10 −5 are considered transients. The system presents two kinds of transients: spirals which lead to X * 1,2 and have negative LEs (black line segments, Figure 1a), and chaotic transients with positive maximal LE (light green line segments, Figure 1a). The spectrum of the local finite-time LEs, Λ = {λ 1 , λ 2 , λ 3 } (Figure 1b), is determined numerically from the system equations with the Wolf algorithm. Considering the lifetime classification mentioned before, irregular or regular trajectories, which continue their dynamics beyond T = 10, 000, are considered in this paper attractors as follows: periodic attractors (see red vertical lines in Figure 1a; stable cycles of the RF system, are not considered in this paper, see [24,25]), quasiperiodic attractors (a quasiperiodic attractor is presented in Figure 2b), hidden chaotic attractors (blue line segments, Figure 1a), and hidden strange nonchaotic attractors (dark green line segments, Figure 1a). Quasiperiodic attractors and irregular transients are merged in the parameter space. There is little distinction between SNAs (in this paper, the attractor H corresponding to b = 0.2875) and chaotic attractors (here, the HCA corresponding to b = 0.28708), as can be seen in the top of Figure 1. They superficially look quite similar, with their structure remaining actually almost unchanged, and even the maximum LEs have the opposite sign. However, as will be discussed next, there are dynamically significant distinctions (see, e.g., [6]).
Equilibria X * 1,2 are stable (stable focus node) for the chosen parameter values, and orbits starting in close neighborhoods will be attracted by X * 1,2 . Equilibria X * 3,4 are unstable (attracting saddles) and orbits starting from close neighborhoods or the attraction basin of X * 3,4 will be attracted by X * 3,4 on the surface determined by the two-dimensional stable manifold; then, after some time, they will exit along the direction of the 1-dimensional unstable manifold.
X * 0 is unstable (repelling focus saddle) for all values of a and b, being globally asymptotically unstable [24], and is a global attractor for the reduced system on x 3 = 0. Therefore, trajectories starting in the plane x 3 = 0 within neighborhoods of X * 0 will diverge to infinity, while trajectories in neighborhoods of X * 0 with x 3 > 0 are attracted by scrolling and by the stable plane of unstable equilibria X * 3 or X * 4 , after which they are repulsed by the unstable directions to infinity. Additionally, depending the initial conditions on the neighborhood of X * 0 , they can be attracted by X * 1,2 along the stable one-dimensional manifold of these equilibria.
Because all equilibria have a pair of complex eigenvalues, all trajectories starting within neighborhoods of equilibria present a scrolling dynamic. Details on equilibria X * 0 and X * 1−4 can be found in [24,25]. Consider next, for a = 0.1, the attractors corresponding to b = 0.98, with the initial condition (0.158, 6 Figure 3f (Appendix B) are presented, respectively. The positiveness of the maximal LE, λ 1 = 0.062, indicates that the attractor is chaotic. This fact is underlined also by the median K ≈ 0.98, close to 1. The broadband of the PSD indicates that the system is chaotic for this value of b. Moreover, the behavior of q as a function of p and the mean-square displacement M are typical for chaotic dynamics (see Appendix A). The RP, with = 0.95 and t ∈ [0, 1000], presents intricate structures that are typical for chaos.
Right after the critical point b = 0.2567, when λ 1 = λ 2 = 0 and λ 3 = −0.314 (Figure 4a), for b = 0.25675 the obtained attractor is quasiperiodic (Figure 2b). For b ∈ [0.2567, 0.2688], the variation in Λ seems to unveil more bifurcations. However, because of the difficulty of analytical study, the bifurcations are not considered here. As for regular motions, K is close to 0, K ≈ −0.001 (Figure 4b). The quasiperiodicity is revealed by the PSD (Figure 4c), whose peaks represent the sub-harmonics (multiples) of a main harmonic generated at b = 0.2567. q as a function of p (Figure 4d) and the evolution of the mean-square displacement M (Figure 4e) show that the attractor is quasiperiodic. Besides the continuous equidistant and parallel off-second diagonal line in the recurrence plot, with = 0.95 and t ∈ [0, 1000] (Figure 4f), there exist another set of interrupted short off-second diagonal segments, parallel with the first set of lines (see the larger circled zoom in Figure 4f). Since all these lines are non-equidistant, the dynamics are quasiperiodic. Note that another set of interrupted off-second diagonal lines seem to be born (see the smaller circled zoom in Figure 4f), suggesting a supplementary frequency. However, this frequency cannot be seen in the PSD (Figure 4c).    [37]. Via the mentioned numerical tools, it is shown that the attractor presents the characteristic of SNAs. The attractor H is hidden because its attraction basin has no connection with neighborhoods of unstable equilibria (see [26] for details) and trajectories starting from small neighborhoods of unstable equilibria X 3,4 and X * 0 are attracted either by the stable equilibria X 1,2 or diverge to infinity, as shown in Figure 2c. For clarity, in Figure 2c only one trajectory from each equilibrium is presented. As can be seen, trajectories starting from a neighborhood of the unstable equilibrium X * 0 tend either to X * 1 (orange plot), to X * 2 (light green plot), or to infinity via the unstable equilibrium X * 3 (black plot) or via the unstable equilibrium X * 4 (grey plot). Trajectories from the neighborhood of X * 3 are attracted either by X * 2 (red plot) or by infinity (brown plot). Trajectories starting from the neighborhood of X * 4 are attracted either by X * 1 (blue plot) or by infinity (brown plot). The attraction basin of H within a parallelipipedic neighborhood of the unstable equilibrium X * 0 is plotted by small red spheres in Figure 5a (note that the system is defined for x 3 > 0 and the axis x 3 is invariant with the reduced equationẋ 3 = −2bx 3 , which has the solution x 3 (t) = e −2bt x 2 (0) [24]). In Figure 5b, one can see a zoomed cubic area around the unstable equilibrium X * 0 , which underlines the fact that, for the considered resolution, there exists an empty cone-like neighborhood with a radius of about 0.05, without initial conditions (red points) of H, indicating that the attraction basin of H has no connection with small neighborhoods of the unstable equilibrium X * 0 . From Figure 5c, one can see that beside the red points (initial conditions) leading to the hidden irregular attractor H, there are two other kinds of initial points which lead to the attractive equilibria X * 1,2 (blue and green plots) [37]. (c) the same zoomed area presenting, as a supplement, the initial conditions leading to the stable equilibria X * 1,2 (blue and red plot).

SNAs of the RF System
Next, it is verified that the hidden attractor H verifies the characteristics of SNAs. For this purpose, the '0-1' test, PSD, and RP are used. In Table 1, the characteristics of the considered attractors of the RF system are presented comparatively.  As shown in Figure 6, the LEs are negative for a small but connected interval of b, V b * (Figure 6a); the oscillations of LEs outside V b * indicate presumably attractors' coexistence; the median K has an intermediate value between 0 and 1, K = 0.62 (Figure 6b); the few peaks in the PSD (F in Figure 6c) not covered by the broadband in the PSD indicate that the attractor is neither quasiperiodic nor chaotic. To note that in [38] it was shown that the PSD of SNAs is discrete, but very dense. Additionally, q as a function of p has a graph between a perfect disc (or circle), for the case of regular behavior, and the Brownian motion (Figure 6d), for chaos; the mean-square displacement M presents an oscillating but bounded increase (Figure 6e); the RP, with = 0.95 and t ∈ [0, 500], presents many non-equidistant off-second diagonal lines parallel with the second diagonal beside several isolated points.
Moreover, the strangeness character, which refers to the nontrivial complicated fractal geometry of the attractor, can be seen from, e.g., the dotted plot in the zoomed area D in Figure 2c, which, besides the non-differentiability of the trajectory, reveals the repeating isolated dust-like points typical for fractals and the correlation dimension D 2 (Appendix C) is D 2 = 1.2.
Since the system is time-continuous, the time series are generally oversampled. This happens because most numerical integrators have variable step sizes, which can generate extremely large numbers of points within some critical trajectory zones. Therefore, before applying the '0-1' test, the underlying time series has to be resampled. Otherwise, it is possible to obtain small (possible zero) values of K, even if the underlying attractor is chaotic [39]. A simple way to resample in Matlab is to use the function resample.
Because of the negativity of the maximal LE on the HSNA H, trajectories do not diverge by separating from each other and eventually coincide, as presented in Figure 7. Comparatively, the first components x 1 (red plot) and y 1 (blue plot) of two trajectories x and y, starting from two nearby initial conditions x 0 and x 0 + (1 × 10 −5 , 1 × 10 −5 , 0) for the SECA corresponding to b = 0.98 (Figure 7a), the quasiperiodic attractor corresponding to b = 0.263 (Figure 7b), and the HSNA H (Figure 7c) are considered. In the case of coupled SNAs, this property characterizes a robust synchronization, in the sense that the non-dependence on initial conditions is expected to synchronize the SNAs quickly, even when they operate with different sets of initial conditions [40].    Figure 7a shows that the difference between the two chaotic trajectories separate-see Figure 7b,c, with details D 2 and D 3 -suggesting that both quasiperiodic trajectories and both trajectories in the case of HSNA present periodic evolution while remaining close to each other. Note that this property of the H is not so evident as for SNAs in quasiperiodically driven discrete systems.
Because the RPs do not show the differences between the considered attractors clearly enough, in Figure 8 the histograms of the considered attractors are considered. As can be seen, the HSNA H can be viewed as being situated between the quasiperiodic attractors and the chaotic attractors. Hidden SNAs, such as H, also exist for b ∈ (0.28743, 0.2876) (see Figure 1).

Conclusions and Discussion
In this paper, it has been shown numerically that the attractor corresponding to b = 0.2875 presents the characteristics of SNAs. For this purpose, the following tools were utilized: Lyapunov exponents, the '0-1' test, power spectral density, recurrence plots, histograms, and the difference between close trajectories. The considered attractor H has all negative LEs, the median K value is between 0 and 1 (K ≈ 0.62), the graph of p and q presents a disk-like structure, the mean-square displacement M has an oscillatory bounded growth, and the recurrence plots present non-equidistant off-second diagonal lines parallel to the second diagonal and also several isolated points. Additionally, contrarily to chaotic attractors, it is verified numerically that the distance between two close trajectories on the HSNA H does not change significantly over time. The fact that LEs look almost similar for the HSNAs and regular attractors could be explained by some inaccuracy in LE calculations for the multi-stable regime observed here. On the other hand, differences between the HSNA and regular attractors of this system are clearly underlined by the 0-1 test, recurrence plot, and power spectral density, which are different.
Even in dynamics do not have explicit quasiperiodic forcing, it is shown numerically that the motion of the RF system could appear on strange nonchaotic attractors. Note that discrete dynamical systems, without a quasiperiodic driving force where the motion is on strange nonchaotic attractors, are considered in [41].
Some open problems deserve further analysis. Thus, the obtained results in this paper regarding the SNA characteristics of H should be considered subject to possible numerical errors due to system's strong dependence on the numerical integrators and time integration interval [0, T], here set to T = 10,000 by using the Matlab options function; strong sensibility on the parameters (a perturbation of order 1 × 10 −4 can move the trajectory to different attractors); initial conditions; etc. Probably due to these reasons, the conclusions regarding the SNA characteristics of H are not so evident, just like the cases of quasiperiodic-driven systems. However, the existence of this hidden strange nonchaotic attractor is verified numerically for a connected parameter interval. Therefore, the attractor H would not be just simply a numerical "accident".
Using different numerical methods to integrate the system, the results could be somewhat different but consistent: there exist small parameter intervals where the RF system has HSNAs.
The usual transition quasiperiodicity → SN A → chaos does not seem to be characteristic for this system. Indeed, Figure 1 shows that after the interval of b, which generates HSNAs, it follows a parameter window where, after chaotic (with maximal positive LE) or nonchaotic (with negative maximal LE) transients, the system enters into a regular dynamics domain (stable equilibria X * 1,2 ). Note that in quasiperiodically forced dynamical systems, the values of the Lyapunov exponents are negative for both quasiperiodic attractors (tori) and SNAs, and in this case the Lyapunov exponents fail to detect transitions from quasiperiodic dynamics to SNAs [42].
On the other hand, an open problem here is the presumably weak chaos. Weak chaos means that even though two trajectories are initially close, they will diverge to be weaker, and their divergence remains finally bounded. In other words, trajectories are not simply regular nor fully chaotic. The separation of nearby trajectories in the case of weak chaos implies that the corresponding maximal Lyapunov exponent is zero [43,44]. Sometimes, this situation appears for some large parameter intervals (see, e.g., the light green vertical bands in Figure 9a). Therefore, the reliability of the '0-1' test [45], used for this system, could be interpreted from the perspective of the weak chaos. (Figure 9), because the maximal LE, λ 1 , is 0 (Figure 9a), the oscillations of K (Figure 9b or Figure 4b) could be due to the existence of weak chaos. For b ∈ [b 4 , b 5 ], K increases but cannot reach the value corresponding to chaos, 1, and the maximal LE is only slightly larger than 0. A clear image for b ∈ [b 5 , 0.288] can be seen in Figure 6a,b. Actually, even the advantage of the '0-1' test is due to the fact that it is binary: 'white' (0) and 'black' (2).
The 'gray'-like situations, such as values 'close' to 1 or 0 for the RF system, require a careful supplementary analysis to clearly state that H is non-chaotic. Additionally, as Figure 7 show, even when two close trajectories remain at a constant distance to each other, the difference between them seems not to vanish at the end of a relatively large integration interval, unlike in the case of SNAs in quasiperiodically driven discrete systems, where, under some circumstances, trajectories eventually coincide, which represents a characteristic of SNAs: robust synchronization.
The use of the recurrence quantification analysis [46] for the detection of structural changes in the dynamics of complex nonlinear systems could give a supplementary clarification of H.
Searching SNAs in other continuous time systems could offer more exciting results. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from authors.
Acknowledgments: Figure 5 was prepared with help from Paul Bourke.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The '0-1' Test
The '0-1' test, which has its roots in [47] and was developed in [48] (see also [39] or [42]), is designed to distinguish chaotic behavior from regular behavior in deterministic systems, which could be used to identify SNAs as well. However, the '0-1' test has been used to identify some SNAs (see, e.g., [42]). The test is easy to implement and does not need system equations, with the input being a time series. Consider a dynamical system with a onedimensional observable data set, constructed from a time series of the underlying system, φ(j), j = 1, 2, ..., N. The '0-1' test is based on a theorem, which states that a nonchaotic motion is bounded, while a chaotic dynamic behaves like a Brownian motion [47]. (1) First, for c ∈ [0, 2π], compute the translation variables p and q by [48]: for n = 1, 2, . . . , N. The choice of c represents an important and sensible algorithm variable (see, e.g., [39], where c is recommended on the interval [π/5, 4π/5]). (2) To determine the growths of p and q, here we used the mean-square displacement M: where n N (in practice, n = N/10 would be a good choice).
The asymptotic growth rate K is defined as follows: Because of the occurrence of resonances for some isolated values of c (where K c could be larger), the median of the computed values of K c (generally, 100 values), denoted by K, is used, since the median is robust against outliers associated with resonances [48]. If the underlying dynamics are regular (i.e., periodic or quasiperiodic), then K ≈ 0; if they are chaotic, then K ≈ 1. For SNAs, K takes an intermediate value between 0 and 1.
In Figure A1a, for the cases of the logistic map x n+1 = rx n (1 − x n ) and the GOPY map [2] x n+1 = 2a tanh(x n ) cos(2πθ n ), θ n+1 = θ n + ω, the plots of q versus p are shown, while in Figure A1b M is presented as a function of n. Figure A1a(i),b(i) represent the regular orbit of the logistic map for r = 3.55; Figure A1a(ii),b(ii) present the chaotic orbit for r = 4, while Figure A1a(iii),b(iii) present the SNA of the GOPY map with a = 1.5 and ω = ( (iii) presents an SNA of the GOPY map x n+1 = 2a tanh(x n ) cos(2πθ n ), θ n+1 = θ n + ω, with a = 1.5, and ω = (

Appendix B. Recurrence Plots
Recurrences in dynamical systems have been analyzed since Poincaré [49], who showed that under certain conditions, orbits of a bounded dynamical system return arbitrarily close to each former points of its own route. However, the time of this return can be arbitrarily long. A tool to visually analyze this recurrence is provided by recurrence plots (RPs), introduced by Eckmann et al. in [50]. Given an autonomous dynamical system and a trajectory {x i } N i , x i ∈ R n , the RP is based on the matrix: Here, is a predefined threshold, · is a norm (Euclidean or maximum norm), and Θ is the Heaviside function. The -recurrent (close) points give "1" while separate points give "0". In a two-dimensional plot, the value is 1 when a point at position j falls into the neighborhood of radius of point at i, and 0 when the point at j does not enter the neighborhood of points at i; they are depicted as black and white dots, respectively. In this way, a visual representation of the system dynamics is provided. Generally, the threshold ε should be chosen small enough. However, in the case of noise influence, a larger threshold is necessary. At some time moment, the state of the system will recur as close as one wishes ( threshold) to a former state. RPs can help find transitions. Diagonal lines correspond to trajectories passing in the same region of the phase space at different times. The RPs of periodic trajectories are equidistant parallel to the second diagonal lines (with slope of 45 • ), indicating some determinism or periodicity (see Figure A2a for the periodic logistic map, with r = 3.4, and threshold = 0.05). The vertical distances between lines in the RP (perpendicular to the horizontal axis) represents the period. The RP of a quasiperiodic trajectory consists of parallel, off-second diagonal lines with different distances between them, reflecting the existence of different return times. If many intricate structures (clusters) disposed on interrupted parallel lines appear, the RP indicates a chaotic trajectory. Now, the distance between the parallel lines is not constant (see Figure A2b for the chaotic logistic map, with r = 4, and threshold = 1.2). For details regarding other RP parameters (embedding dimension, time delay, etc.), see [51][52][53]. For the use of the RP to study SNAs in the forced logistic map, see [54].

Appendix C. Correlation Dimension
Let an d-dimensional embedding space, x 1 , x 2 , ..., x m a time series and y 1 , y 2 , . . . , y N , with N = m − d + 1, and the delay coordinates y j = (x j , . . . , x j+d−1 ), j = 1, . . . , N. The correlation integral for a given r > 0 is defined as follows: where H is the Heaviside step function and |.| is the Euclidean distance. Here, the H counts the number of pairs (x i , x j ), satisfying the condition |x i − x j | < r.
If C(r) scales like C(r) = r D 2 , then: is called the correlation dimension, the straight line of log-log curve of r and C(r). The G-P algorithm, created by Grassberger and Procaccia [55], is one of the most used algorithms to determine D 2 numerically, and uses the embedding theory together with the idea of phase space reconstruction.