Dynamical robustness of discrete conservative systems: Harper and generalized standard maps

The dynamics of the ergodic or mixing systems can be explained via Boltzmann–Gibbs (BG) statistical mechanics. For these systems, exponential and Gaussian distributions are appropriate forms to describe the limit probability behavior of the relevant variables of the system under consideration. Since these distributions maximize the BG entropy, the reason behind the occurrence of such distributions can be explained by the central limit theorem (CLT). On the other hand, due to the ergodicity breakdown and strong correlations among the random variables observed for some systems for some intervals of system parameters, BG statistical mechanical approaches fail to describe the dynamics of those cases and the CLT is not valid anymore. It has been shown in recent years that a generalization of the CLT is possible and the limit probability distributions seem to converge to a q-Gaussian distribution [1–6] for a class of systems with certain correlations. Like the role of the Gaussians in BG statistics, q-Gaussian distributions maximize the nonadditive entropy (${S}_{q}=k\left(1-{\sum }_{i}{p}_{i}^{q}\right)/\left(q-1\right)$) and constitute the basis of nonextensive statistical mechanics [7, 8]. The nonextensive statistical mechanical framework provides a more general picture by recovering the BG statistical framework as a special case (q → 1) and this generalization seems to be the appropriate tool for explaining the statistical mechanical behavior of a wide range of systems where BG statistics is known to fail.

In recent years, domains of validity of these two statistical mechanical frameworks have been shown by utilizing the rich phase space behaviors exhibited in area-preserving maps [9–11] (namely, the standard map and the web map). As the amount of nonintegrability increases with the increment of the map parameter for these Hamiltonian systems, chaotic and regular behavior regions may coexist in the phase spaces of these maps for specific parameter values. In the chaotic regions, the system is ergodic and iterates of the chaotic trajectories display uncorrelated behavior while wandering throughout the allowed region in the phase space. In this case, the limit probability distribution of the sums of iterates of the map converges to a Gaussian consistently with assertions of the BG statistical mechanical framework. On the other hand, iterates of trajectories starting from the initial conditions located inside the nonergodic stability islands exhibit a strongly correlated behavior. For these initial conditions, it has been shown for the standard map [9, 10] and the web map [11] that the limit probability distribution converges to a q-Gaussian with a specific q = 1.935 ± 0.005 value. When the entire system is modeled by using the initial conditions coming from both chaotic sea and stability islands for some specific parameter values, the limit probability distribution is obtained as a linear combination of a Gaussian and a q-Gaussian with q = 1.935 ± 0.005; the q-Gaussian distribution maintains its existence together with a Gaussian even for a large number of iteration steps. This limit distribution seems very robust from a statistical mechanical point of view since the same distribution appears for the initial conditions selected from the stability islands of different maps. Any novel understanding of such systems would no doubt be important if we consider the role of the area-preserving maps in physics and in the development of the chaos theory [12]. Many physical systems such as magnetic traps [13], electron magnetotransport in classical and quantum wells [14], and particle accelerators [15] can be modeled by using the standard map as a first approximation. In addition, a linear combination of the standard map and the web map can be used for modeling many physical systems [16, 17]. As the area-preserving maps are of great importance in different branches, whether the common limit probability behavior observed for the initial conditions of the stability islands of the standard map and the web map is a universal behavior for all area-preserving maps is a very intriguing research question that deserves to be investigated. In this study, to test the robustness of the q-Gaussian distribution with q = 1.935 ± 0.005, we analyze the limit behavior of the sums of the iterates of the Harper map [18], which models transport phenomena in deterministic chaotic Hamiltonian systems. In addition to the Harper map, we also define a new generalized form of the standard map, namely a z-generalized standard map, to create independent and unique area-preserving maps exhibiting different phase space dynamics.

In order to statistically characterize these area-preserving maps we choose various map parameter values for each system where the phase spaces display different behaviors. As the phase spaces of these scenarios provide rich observations by exhibiting chaotic and regular behaviors at the same time for specific parameter values, we numerically investigate the limit probability distribution of the entire system and visualize how the limit distributions vary according to the map parameter value. For each case we obtain phase space portraits given in this paper by iterating 40–50 randomly chosen initial conditions T = 5 × 10³ times. In order to quantify the trajectory behaviors seen in the phase portraits we calculate the finite-time largest Lyapunov exponent, FLLE (λ), as discussed in [19] by using the Benettin algorithm [20] for each initial condition randomly chosen from the entire phase space. Lyapunov exponents are calculated over T = 5 × 10⁵ time steps using M = 5 × 10⁵ initial conditions and the Lyapunov spectra of scenarios are portrayed via color maps in order to reveal the regions with different behavior. Since the chaotic regions and the stability islands exhibit largely positive and nearly zero FLLE values, respectively, this separation of the phase space regions allows us to distinguish the portions of the phase space where the system appears to be ergodic and nonergodic [9]. Chaotic trajectories presenting largely positive FLLE values diverge exponentially in the allowed region of the phase space and these trajectories spread into this region with apparently random behavior. For the strongly chaotic regions, the system exhibits mixing property and ergodic behavior. On the other hand, trajectories located inside the stability islands, which can be periodic or quasi-periodic, present nearly zero FLLE values (λ ≈ 0) and the system is nonergodic in these regions. Although the argument for the ergodicity of the chaotic trajectories has been verified for numerous dissipative [23, 24] and conservative [9–11] systems, we came across a contrary situation for sticky chaotic regions that may arise in several z-generalized standard map systems by exhibiting nonergodic behavior for finite observation times. This unexpected behavior will be further discussed later on when analyzing the z-generalized standard map.

Since the ergodic and nonergodic portions of the phase space require different statistical mechanical approaches, we can investigate the limit probability distributions of the variables of the systems to determine the domains of validity of the BG and of the nonextensive statistical mechanics. In the spirit of the central limit theorem, for the limit probability distribution characterization, we define the variable

$$\begin{equation}y=\frac{1}{\sqrt{T}}\sum _{i=1}^{T}\left({x}_{i}-\langle x\rangle \right)\end{equation} \tag{ 1 }$$

where x is the variable of the map, ⟨⋯⟩ denotes averaging over a large number of iterations T and a large number of randomly chosen initial conditions M, i.e., $\langle x\rangle =\frac{1}{M}\frac{1}{T}{\sum }_{j=1}^{M}{\sum }_{i=1}^{T}{x}_{i}^{\left(j\right)}$. It has been previously shown that for arbitrary values of the parameter K of the standard map [9, 10], that the probability distribution of these sums (equation (1)) can be modeled as

$$\begin{equation}{P}_{q}\left(y;{\mu }_{q},{\sigma }_{q}\right)={A}_{q}\sqrt{{B}_{q}}{\left[1-\left(1-q\right){B}_{q}{\left(y-{\mu }_{q}\right)}^{2}\right]}^{\frac{1}{1-q}},\end{equation} \tag{ 2 }$$

which represents the probability density for the initial conditions inside the vanishing Lyapunov region (q ≠ 1), where μ_q is the q-mean value, σ_q is the q-variance, A_q is the normalization factor, and B_q is a parameter that characterizes the width of the distribution [25]:

$$\begin{equation}{A}_{q}=\begin{cases}\frac{{\Gamma}\left[\frac{5-3q}{2\left(1-q\right)}\right]}{{\Gamma}\left[\frac{2-q}{1-q}\right]}\sqrt{\frac{1-q}{\pi }},\hfill & q{< }1\hfill \\ \frac{1}{\sqrt{\pi }},\hfill & q=1\hfill \\ \frac{{\Gamma}\left[\frac{1}{q-1}\right]}{{\Gamma}\left[\frac{3-q}{2\left(q-1\right)}\right]}\sqrt{\frac{q-1}{\pi }},\hfill & 1{< }q{< }3\hfill \end{cases}\end{equation} \tag{ 3 }$$

$$\begin{equation}{B}_{q}={\left[\left(3-q\right){\sigma }_{q}^{2}\right]}^{-1}.\end{equation} \tag{ 4 }$$

The q-mean value and q-variance are defined by (see [25] for the continuous version):

$$\begin{equation}{\mu }_{q}=\frac{\sum _{i=1}^{N}{y}_{i}{\left[{P}_{q}\left({y}_{i}\right)\right]}^{q}}{\sum _{i=1}^{N}{\left[{P}_{q}\left({y}_{i}\right)\right]}^{q}},\end{equation} \tag{ 5 }$$

$$\begin{equation}{\sigma }_{q}^{2}=\frac{\sum _{i=1}^{N}{y}_{i}^{2}{\left[{P}_{q}\left({y}_{i}\right)\right]}^{q}}{\sum _{i=1}^{N}{\left[{P}_{q}\left({y}_{i}\right)\right]}^{q}},\end{equation} \tag{ 6 }$$

though we have considered these variables as fitting parameters.

The q → 1 limit recovers the Gaussian distribution ${P}_{1}\left(y;{\mu }_{1},{\sigma }_{1}\right)=\frac{1}{{\sigma }_{1}\sqrt{2\pi }}\mathrm{exp}\left[-\frac{1}{2}{\left(\frac{y-{\mu }_{1}}{{\sigma }_{1}}\right)}^{2}\right]$. In the analyses of the limit probability distributions of such maps with various values for the map parameters, we randomly chose a large number of initial conditions, larger than 3 × 10⁷, from the entire phase space, and numerically calculate the limit distribution of equation (1) using T = 2²² iteration steps in order to obtain a good statistical description of the systems. These values are determined in accordance with recent works [9, 10] where they have been shown to be optimal by considering the required computational times and convergence of the obtained probability distributions.

This paper is organized as follows: firstly, the Harper map and the z-generalized standard map are introduced in sections 2 and 3; then the results of the numerical calculations are discussed in section 4. Finally, we conclude in the last section.

In order to investigate transport phenomena in deterministic chaotic Hamiltonian dynamics, a two-dimensional area-preserving map has been proposed in [18] using a time-dependent Hamiltonian of the form

$$\begin{equation}H\left(x,p,t\right)=-{V}_{2}\mathrm{cos}\left(2\pi p\right)-{V}_{1}\mathrm{cos}\left(2\pi x\right)\tau \sum _{n=-\infty }^{\infty }\delta \left(t-n\tau \right).\end{equation} \tag{ 7 }$$

This is called the kicked Harper model and has several applications in physics [26–28]. If one integrates equation (7) over one period τ of the kicking potential, the kicked Harper map can be obtained easily as follows:

$$\begin{equation}\begin{aligned}\hfill {p}_{i+1}& ={p}_{i}-{\gamma }_{1}\mathrm{sin}\left(2\pi {x}_{i}\right)\hfill \\ \hfill {x}_{i+1}& ={x}_{i}+{\gamma }_{2}\mathrm{sin}\left(2\pi {p}_{i+1}\right)\hfill \end{aligned}\end{equation} \tag{ 8 }$$

where p and x are taken as modulo 1 and γ_j = 2πV_jτ. In the present paper, we only consider the special case γ₁ = γ₂ ≡ γ.

The Hamiltonian of the standard map is given by [29]

$$\begin{equation}H=\frac{1}{2}{p}^{2}-K\mathrm{cos}\left(x\right)\sum _{n=-\infty }^{\infty }\delta \left(\frac{t}{T}-n\right)\end{equation} \tag{ 9 }$$

where p and x are the momentum and position of a particle, respectively, and the periodic sequence of δ-pulses has the period T = 2π/ν. The equations of motion derived from this Hamiltonian enables us to write down the momentum and position variables at the nth and (n + 1)th kicks, from where the original standard map is obtained [30]. One possible way of generalizing the original standard map is to modify the kicked term as (K/z)cos(zx), which results in defining the z-generalized standard map as follows:

$$\begin{equation}\begin{aligned}\hfill {p}_{i+1}& ={p}_{i}-K\mathrm{sin}\left(z{x}_{i}\right)\hfill \\ \hfill {x}_{i+1}& ={x}_{i}+{p}_{i+1}\hfill \end{aligned}\end{equation} \tag{ 10 }$$

where p and x are taken as modulo 2π, K is the map parameter that controls the amount of nonintegrability of the system, and z is an integer; z = 1 recovers the usual standard map. Since K = 0 case is an integrable, there is no arbitrary z value that can make the system completely chaotic. However, in principle, for a small deviation from the K = 0 case, there must be a large value of z which makes the system chaotic. With different z values we define unique systems exhibiting specific phase-space dynamics. In this paper we investigate the phase-space behaviors and the limit distributions for various z values with K = 0.2 and K = 0.6 parameters by using the numerical calculations introduced in the introduction section. In order to present a clear evolution of the phase portraits and the limit distributions according to the z generalization term, we analyze z = 3, z = 5, z = 40 generalized systems for K = 0.2, and z = 3, z = 4, z = 15 systems for K = 0.6. For both parameter values we also analyze the z = 1 system (the original standard map) in order to study how the z-generalization modifies the dynamical and statistical behavior of the system.

As discussed in [18], in the Harper map phase-space plots, a separatrix defined by H₀(p, x) = 0 appears and forms a square symmetry. When γ > 0, this separatrix is destroyed and becomes a mesh of finite thickness inside which the dynamics is chaotic. This region grows as γ increases, making the region of the regular motion shrink. This behavior is evident from the first column of figure 1 for some representative γ values. Then, in the second column, one can see the Lyapunov diagram for the same values of γ. The genesis and increasing domination of the chaotic sea can be clearly seen as γ increases. Finally, the last column exhibits the corresponding probability distributions. As discussed extensively in [10, 11], these distributions can be well approximated by a linear combination of one Gaussian and one q-Gaussian distribution, namely,

$$\begin{equation}P\left(y\right)={\alpha }_{{q}_{1}}{P}_{{q}_{1}}\left(y;{\mu }_{{q}_{1}},{\sigma }_{{q}_{1}}\right)+{\alpha }_{{q}_{2}}{P}_{{q}_{2}}\left(y;{\mu }_{{q}_{2}},{\sigma }_{{q}_{2}}\right)\end{equation} \tag{ 11 }$$

where q₁ = 1.935 ± 0.005 and q₂ = 1. In equation (11), the contribution ratios ${\alpha }_{{q}_{1}}$ and ${\alpha }_{{q}_{2}}$ can be evaluated from the phase-space occupation ratios of the initial conditions located in the stability islands and chaotic sea detected from the Lyapunov color map, respectively. Therefore, these parameters are not fitting parameters, but determined directly from the dynamics of the system. The q-Gaussian distribution with q = 1.935 ± 0.005 originates from the initial conditions of the stability islands, and the initial conditions from the chaotic sea contribute to the Gaussian distribution. The obtained results for the parameters are given in table 1. These results clearly show that the Harper map is in the same universality class of the standard map and the web map, and therefore provides another argument pointing to the robustness of the q-Gaussian with q = 1.935 ± 0.005.

Zoom In Zoom Out Reset image size

Table 1. Obtained results for the Harper map for four representative values of γ.

	γ = 0.05	γ = 0.12	γ = 0.15	γ = 0.9
q1	1.935	1.935	1.935	1.935
q2	1	1	1	1
${\alpha }_{{q}_{1}}$	1	0.9960	0.9713	0
${\alpha }_{{q}_{2}}$	0	0.0040	0.0287	1
${B}_{{q}_{1}}$	65 × 105	25 × 105	25 × 105	—
${B}_{{q}_{2}}$	—	5.3	5.7	7.0

Now, we can concentrate on the z-generalized standard map. Phase-space portraits of representative cases, their corresponding Lyapunov diagrams, and the limit probability distributions are given in figure 2 for K = 0.2 and in figure 3 for K = 0.6. They enable the visualization of how the dynamics of these systems change according to the increment of the z term. Surprisingly, we noticed that these distributions cannot be modeled using equation (11). Instead, we verified that, for this system, the obtained distributions happen to be well approximated by a linear combination of three q-Gaussian distributions, namely,

$$\begin{equation}P\left(y\right)={\alpha }_{{q}_{1}}{P}_{{q}_{1}}\left(y;{\mu }_{{q}_{1}},{\sigma }_{{q}_{1}}\right)+{\alpha }_{{q}_{2}}{P}_{{q}_{2}}\left(y;{\mu }_{{q}_{2}},{\sigma }_{{q}_{2}}\right)+{\alpha }_{{q}_{3}}{P}_{{q}_{3}}\left(y;{\mu }_{{q}_{3}},{\sigma }_{{q}_{3}}\right).\end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The probability distributions given in figures 2 and 3 can be well approximated by this linear combination. The relevant parameter values are given in table 2.

Table 2. Obtained results for the z-generalized standard map for representative values of z and K.

	K = 0.2		K = 0.6
	z = 3	z = 5	z = 3	z = 4
q1	1.935	1.935	1.935	1.935
q2	1.40	1.55	1.40	1.55
q3	1	1.45	1	1
${\alpha }_{{q}_{1}}$	0.963	0.515	0.265	0.114
${\alpha }_{{q}_{2}}$	0.025	0.340	0.565	0.300
${\alpha }_{{q}_{3}}$	0.012	0.145	0.170	0.586
${B}_{{q}_{1}}$	60 000	25 800	25 300	2000
${B}_{{q}_{2}}$	0.389	0.072	0.0215	0.0039
${B}_{{q}_{3}}$	0.000 15	0.004	0.0455	0.0025

It can be seen in figure 2 that the phase spaces of the z = 1 and z = 40 systems are entirely occupied by nonergodic stability islands and ergodic chaotic sea, respectively. Conformably with these phase space behaviors, the probability distribution of equation (1) obtained for the initial conditions randomly chosen from the entire phase space is well fitted by a Gaussian when the system is ergodic-like and by a q-Gaussian with q ≃ 1.935 when the system is nonergodic-like. In figure 3, for K = 0.6, the occurrence of a Gaussian as a limit distribution of the initial conditions chosen from the ergodic phase space is also verified for the z = 15 system whose phase space is occupied by the chaotic sea. For the K = 0.6 parameter value of the z = 1 system which corresponds to the original standard map, the phase space consists of both stability islands and the chaotic sea. In accordance with recent works [9, 10], the limit probability distribution of equation (1), given in figure 3, is obtained as a linear combination of a Gaussian arises from the initial conditions located in the chaotic sea and a q-Gaussian with q = 1.935 ± 0.005 arises from the initial conditions located in the stability islands. By considering the regions of different behaviors that the phase space consists of, these limit distributions are of course expected due to the ergodic/nonergodic behavior of the related phase space. Here contribution ratios of each term in the linear combination are detected using the Lyapunov spectrum, and therefore they are not fitting parameters but determined directly from the dynamics of the system.

We come across with a more complicated limit behavior for the other parameter values of the generalized systems with the K = 0.2 and K = 0.6 cases. Even though the phase spaces of these systems, given in figures 2 and 3, consist of both stability islands and the chaotic sea similar to the original standard case, the obtained limit distributions exhibit a three-component behavior which can be modeled by equation (12). For each system mentioned above, the relevant parameter values of the probability distributions, shown in figures 2 and 3, are given in table 2. We see that an unexpected third component is obtained as a q-Gaussian with different q values for (K = 0.2, z = 3), (K = 0.6, z = 3), and (K = 0.6, z = 4) cases. More surprisingly, the obtained probability distribution of the (K = 0.2, z = 5) case consists of three q-Gaussians. In order to explain these interesting observations, we have to focus on the requirements for the occurrence of the q-Gaussians and the phase space behavior of the z-generalized standard map.

With the present z-generalization, we actually create systems with different phase space behaviors. As it can be seen from the phase spaces given in figures 2 and 3, for typical values of K, the increment of the z term increases the amount of nonintegrability of the system and chaotic behavior may occur for smaller K values compared to the original standard map case. With the increasing nonintegrability [12], the stability islands which are actually survived Kolmogorov–Arnold–Moser (KAM) tori dissolve according to the KAM theorem [31] and the Poincaré–Birkhoff theorem [32].

According to the Poincaré–Birkhoff theorem, as a result of the resonances arising with increased nonintegrability, each torus dissolves into an alternating sequence of a series of hyperbolic–elliptic points depending on its winding number. Elliptic orbits occur around each elliptic point and in-sets and out-sets of hyperbolic points surround these elliptic orbits. Each elliptic–hyperbolic points series and in-sets and out-sets of hyperbolic points constitute a resonance. As the nonintegrability increases, resonances begin to overlap and destroy surviving tori that were in the region between them [12]. Chaotic behavior occurs with complex tangle structures created by in-sets and out-sets of the hyperbolic points. Homoclinic and heteroclinic tangles surround the elliptic orbits without intersecting them and a chaotic trajectory spread throughout this tangle structure by exhibiting mixing behavior [33]. Dissolution of tori and enlarging chaotic behavior create archipelagos centered around the elliptic points in the phase space. The outer stability islands of the archipelagos continue to dissolve according to the KAM theorem and the Poincaré–Birkhoff theorem and chaotic behavior occupies a larger portion of the phase space. It is important to note that this formation can be used to explain the occurrence of chaotic behavior in Hamiltonian systems and the chaotic sea observed in the phase space corresponds to a single constant energy region of the Hamiltonian. By considering the resonance structure discussed above, we can conclude that complex homoclinic and heteroclinic tangles should stick around the elliptic orbits. This sticky behavior exists in the original standard map but chaotic trajectories spend less time in these sticky regions before escaping into the chaotic sea where they can wander throughout apparently randomly. When we modify the sine term in standard map as given in equation (10), we radically change the resonance behaviors and the alternating sequence of the series of hyperbolic–elliptic points observed in the original case. As it can be seen from archipelagos in figures 2 and 3, each generalized system has different elliptic–hyperbolic points organization with periodicity related to the z term in the sine function. Starting from the integrable case, the case that is common for all z values, tori located in resonance structures dissolve because of the overlaps. The sine term of the generalized standard map affects the number of hyperbolic–elliptic points and their positioning in the phase space and as a result of this effect resonances become stronger with an increasing z term for a constant K value. Thus, more in-sets and out-sets of hyperbolic points create more complex tangle structures around stability islands. As these tangles tend to surround and stick around the stability islands, with more complicated structure chaotic trajectories spend more time covering these tangle structures. Even though these sticky regions are connected to a strongly chaotic sea, as can be seen clearly from figures 2 and 3, an initial condition starting from the sticky region may spend most of its time in that region and may escape into the strongly chaotic sea after unpredictable time steps. The statistical effectiveness of the sticky regions is presented in the probability distributions obtained for the z > 1 systems, i.e. especially for the (K = 0.2, z = 5) case which is explained in detail below. A chaotic trajectory should visit both the chaotic sea and each sticky region, because the chaotic sea we see in the phase space is actually an allowed energy surface of the original Hamiltonian system which is created by the dissolution of the constant energy surface tori. As a chaotic trajectory wanders throughout the allowed energy region apparently randomly, its behavior displays a mixing property and the system is said to be ergodic in that region.

When the phase space is fully occupied by the strongly chaotic sea, the whole system is ergodic and Gaussian distributions are obtained, which indicates the validity of the BG statistical framework. On the contrary, as a consequence of the ergodicity breakdown, q-Gaussians are appropriate distributions that describe the whole system with the phase space entirely occupied by stability islands. As these inferences are verified for the (K = 0.2, z = 1), (K = 0.2, z = 40), and (K = 0.6, z = 15) systems, two-component linear combinations of probability distributions of the (K = 0.6, z = 1) system shows that q-Gaussians arise from the initial conditions selected from the stability islands whereas the Gaussian contribution comes from chaotic trajectories. It is important to note here that ergodicity breakdown alone is not sufficient for the occurrence of the q-Gaussians; indeed, a special type of correlation among random variables is also needed. These requirements are fulfilled for stability islands of area-preserving maps [9–11] and, for chaotic bands, of a band-splitting structure that approaches the chaos threshold of the dissipative logistic map by means of a Huberman–Rudnick-like scaling law [21–23].

We can see from table 2 that the same q-Gaussian contributions with q = 1.935 ± 0.005 come for each system. By taking into consideration that the same distribution is obtained for the initial conditions selected from the stability islands of the original standard map (z = 1) with different parameter values [9, 10], we distinguish stability islands from the phase space by using the Lyapunov diagrams given in figures 2 and 3. Moreover, we verify that the phase-space occupation ratios of the stability islands of each case are exactly the ${\alpha }_{{q}_{1}}$ values given in table 2. If the initial conditions from the stability islands are discarded, we are left with the chaotic sea in the phase space which should be related to the remaining components of the probability distribution given in equation (12). Let us start by analyzing the (K = 0.2, z = 3), (K = 0.6, z = 3), and (K = 0.6, z = 4) cases whose chaotic seas give rise to both a Gaussian and a q-Gaussian. When we look at the chaotic seas in the phase spaces of these systems in more detail, we see that strongly sticky regions which cannot occur in the original standard map arise for all these systems. Due to the resonance behaviors mentioned above, these sticky regions surround archipelagos and create complicated structures which are also connected to the chaotic sea. An initial condition located inside one of these sticky regions may give rise to a chaotic trajectory that stays inside this region for many iteration steps. An initial condition inside the sticky region evolves to cover this region, and after unpredictable iteration steps, this trajectory may escape to the chaotic sea and may wander throughout the sea apparently randomly. In the probability distribution analyses, we use large but finite iteration steps and it seems that during these iteration steps chaotic trajectories located inside the sticky regions cannot display a strong mixing behavior by confining inside these regions and not visiting a large portion of the allowed energy region. Chaotic trajectories that do not enter sticky regions for large iteration steps can freely wander in the large portion of the allowed energy region and give rise to a Gaussian distribution. In principle, if it was possible to leave the system to evolve infinitely, the entire equal energy region would be covered by a single chaotic trajectory. In our observation interval, as some of the chaotic trajectories wander freely, some of them instead spend most of their time in the sticky region. Most probably, this observation may be a plausible explanation of the obtained limit probability distributions of the systems that we investigate here. Observations made for the (K = 0.2, z = 5) case present a more complicated scenario. When we look at the Lyapunov spectrum of this case in figure 2, we see horizontal band-like structures. Based on our observations, the chaotic sea of each band acts as a sticky region in the phase space and all these regions are connected by not having any stability island as a barrier between them. When we analyze the phase-space behavior of the chaotic trajectories by selecting the initial conditions inside the chaotic seas and letting them evolve, we observe that chaotic trajectories do not spread into the allowed energy region randomly as expected from the regular chaotic trajectory like we see in the (K = 0.2, z = 40) and (K = 0.6, z = 15) systems. Instead of spreading into the allowed energy region randomly, chaotic trajectories first cover their bands' chaotic sea and then move into another band. In our observations we see that even for a large number of iteration steps, i.e. T = 10¹⁰, the entire chaotic sea cannot be covered by a trajectory starting from an initial condition. As chaotic trajectories move in the phase space covering firstly one band and then the others, respectively, the mixing property of this system is completely different from that of the original system. This difference seems to be related to the occurrence of two q-Gaussian contributions detected for the limit probability distribution from the chaotic sea.

In order to improve our understanding of the emergence of the q-Gaussians, we have also analyzed the auto-correlation function r_κ defined as follows:

$$\begin{equation}{r}_{\kappa }=\sum _{i=1}^{T-\kappa }\frac{\left({y}_{i}-\langle y\rangle \right)\left({y}_{i+\kappa }-\langle y\rangle \right)}{\sum _{i=1}^{T}{\left({y}_{i}-\langle y\rangle \right)}^{2}}\end{equation} \tag{ 13 }$$

where κ is the time lag, T is the number of iteration steps which constitute the trajectory, and $\langle y\rangle ={T}^{-1}{\sum }_{i=1}^{T}{y}_{i}$ [34]. Iterates of a trajectory are correlated for r_κ≠0 and not correlated for r_κ = 0. For each system given in table 2, we randomly chose a large number of initial conditions from the entire phase space and let the system evolve along T = 2²² iteration steps, which coincides with the iteration number used in the probability distribution computations, starting from these initial conditions. In the computations of r_κ, for each trajectory, we use κ = 10⁵ as a maximum time lag by considering large computation times required for larger time lag values. The obtained results show that all systems exhibit three different tendencies for the auto-correlation function compatible with the probability distributions in the form of equation (12). In figures 4 and 5 we demonstrate the auto-correlation functions of the (K = 0.2, z = 5), and (K = 0.6, z = 4) systems as a function of the logarithm of the time lag respectively to corroborate the explanations given for the nonergodic and nonmixing behavior of the chaotic trajectories. The logarithm of the time lag is used and κ is cut at 10⁵ in order to obtain a better visualization of the oscillatory behavior of the auto-correlation functions. In both figures a common green color is used to indicate the auto-correlation functions of iterates starting from initial conditions located inside one of the stability islands, e.g. (x = 5.042 890 762..., p = 0.154 394 798...) for the (K = 0.2, z = 5) system and (x = 1.617 417 044..., p = 4.837 677 782...) for the (K = 0.6, z = 4) system. As can be seen from the figures, a green curve oscillates around zero with a very large amplitudes and this behavior indicates that the iterates in the stability islands are strongly correlated. In figure 4 we see that red and black curves also oscillate around zero with different amplitudes that are smaller than the amplitude of the green one. Both of these functions are obtained for trajectories whose initial conditions are selected from the chaotic sea of the phase space, i.e. (x = 3.190 574 110......, p = 7.259 267 286... × 10⁻²...) for red and (x = 3.096 157 726..., p = 9.095 801 865... × 10⁻²...) for black. These auto-correlation function behaviors are the three main types of correlations that are observed in the analyses of the (K = 0.2, z = 5) system and through these observations we can say that the chaotic trajectories in the phase space exhibit a correlated behavior which is weaker than the correlation among the iterates of a trajectory in the stability islands. These three types of correlations observed in the phase space together with the nonergodicity of the present system mentioned above fulfill the requirements of the occurrence of the q-Gaussians and are thought to explain the obtained limit probability distribution which is a linear combination of three q-Gaussians. In figure 5 red and black curves are obtained for the initial conditions chosen from the chaotic region in the phase space. As the initial condition (x = 5.654 566 893..., p = 1.627 640 289...) a red curve is located in the sticky region, and as the initial condition (x = 5.668 539 896..., p = 4.509 105 458...) a black curve is located in the strongly chaotic sea. When we look at the figure, we see that the black curve decreases to zero after a short time lag and oscillates around zero by indicating the uncorrelated nature of the iterates of the trajectory. This auto-correlation function behavior is similar to the auto-correlation function obtained for the white noise which was recently shown in reference [24]. On the contrary, the red curve oscillates around zero with large amplitudes like the previous scenario. From these observations we can deduce that chaotic trajectories located inside the sticky regions may show correlations and chaotic trajectories that do not enter into the sticky regions for a long period of iteration steps display the expected uncorrelated behavior of the chaotic trajectories. Also, we obtain the same auto-correlation function behaviors for the (K = 0.2, z = 3) and (K = 0.6, z = 3) systems that exhibit a similar limit probability distribution like the (K = 0.6, z = 4) system. These observations seem to provide an adequate explanation for the obtained limit probability distributions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, our results on the area-preserving maps can be summarized by classifying them into two groups: observations for the stability islands and for the chaotic trajectories. For the Harper map and several z-generalized standard map systems, for a large number of iteration steps, the limit probability distributions coming from the sum of the iterates when the system is started from initial conditions located inside the stability islands seem to converge to a q-Gaussian with a q = 1.935 ± 0.005 value. Each different z value corresponds to a new area-preserving system. Taking into account all of the maps analyzed in this paper and also results of the recent paper on the statistical characterization of the area-preserving web map [11], the main goal of this manuscript is to verify numerically that the limit probability distribution obtained when the system is initiated from initial conditions located in the stability islands is always well approximated by a q-Gaussian with a q = 1.935 ± 0.005 value. Regardless of the magnitude of the phase-space occupation ratios of the stability islands (these ratios are not fitting parameters since they come directly from the Lyapunov spectrum of the system) for various map parameter values, a q-Gaussian with q = 1.935 ± 0.005 maintains its presence together with other distributions and this fact indicates that a q-Gaussian with q = 1.935 ± 0.005 value is a robust limit behavior for the stability islands of the area-preserving maps.

Although the stability islands of the area-preserving maps exhibit the same limit behavior, unexpected observations are made for chaotic trajectories of the different maps. Considering the definitional properties of the chaotic trajectories, e.g. the apparently random behavior and the exponential divergence of initially nearby trajectories, one can suggest that the chaotic trajectories wander freely through the allowed energy region and they spread into this region by displaying mixing property. Under normal circumstances, a single chaotic trajectory rapidly spreads into an allowed region and outlines this region after a few iteration steps. This common behavior of the chaotic trajectories is observed for all the control parameter values of the Harper map. For all cases of the Harper map, when the chaotic region develops, we obtain Gaussian distributions for the limit behavior of the sums of the iterates for the initial conditions started from the chaotic sea, as expected. However, as shown here, some chaotic trajectories cannot exhibit the regular behavior of the chaotic trajectories due to the sticky regions occur around the stability islands for some area-preserving maps. Such a chaotic trajectory may not visit most of the allowed energy region during the observation time and therefore it cannot behave similarly as the chaotic trajectories which do not visit sticky regions during the same period of time. Even though chaotic trajectories exponentially diverge while covering the sticky region, the magnitude of their divergence is much smaller compared to the divergence in the chaotic sea as seen from the Lyapunov spectra. A second q-Gaussian is obtained in the limit probability distributions of systems that exhibit sticky behavior in their phase spaces and this might be explained due to different mixing properties compared to the standard case and the correlated nature of the chaotic trajectories of sticky region. Even though the contribution ratios of these second q-Gaussians cannot be determined directly from the Lyapunov spectra as we did before, they are thought to be as robust as other distributions that make contributions to the limit behavior for long but finite time intervals.

We acknowledge fruitful remarks from two anonymous referees. This work has been partially supported by TUBITAK (Turkish Agency) under the Research Project number 115F492, and by CNPq and Faperj (Brazilian Agencies).