domain: phys nlin

text:
Henok Tenaw Moges

Thanos Manos

Charalampos Skokos

Anomalous diffusion in single and coupled standard maps with extensive chaotic phase spaces

Keywords: dynamical systems, standard map, anomalous diffusion, accelerator modes, generalized alignment index (GALI), Lyapunov exponent

We investigate the long-term diffusion transport and chaos properties of single and coupled standard maps. We consider model parameters that are known to induce anomalous diffusion in the maps' phase spaces, as opposed to normal diffusion which is associated with Gaussian distribution properties of the kinematic variables. This type of transport originates in the presence of the so-called accelerator modes, i.e. non-chaotic initial conditions which exhibit ballistic transport, which also affect the dynamics in their vicinity. We first systematically study the dynamics of single standard maps, investigating the impact of different ensembles of initial conditions on their behavior and asymptotic diffusion rates, as well as on the respective time-scales needed to acquire these rates. We consider sets of initial conditions in chaotic regions enclosing accelerator modes, which are not bounded by invariant tori. These types of chaotic initial conditions typically lead to normal diffusion transport. We then setup different arrangements of coupled standard maps and investigate their global diffusion properties and chaotic dynamics. Although individual maps bear accelerator modes causing anomalous transport, the global diffusion behavior of the coupled system turns out to depend on the specific configuration of the imposed coupling. Estimating the average diffusion properties for ensembles of initial conditions, as well as measuring the strength of chaos through computations of appropriate indicators, we find conditions and systems' arrangements which systematically favor the suppression of anomalous transport and long-term convergence to normal diffusion rates.

Introduction

The investigation of diffusion and transport phenomena in conservative Hamiltonian systems and area-preserving symplectic maps is a topic of intense scientific research for more than 50 years now (see e.g. [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF][START_REF] Rechester | Calculation of turbulent diffusion for the Chirikov-Taylor model[END_REF][START_REF] Rechester | Fourier-space paths applied to the calculation of diffusion for the Chirikov-Taylor model[END_REF][START_REF] Cary | Statistical characterization of periodic, area-preserving mappings[END_REF][START_REF] Meiss | Correlations of periodic, area-preserving maps[END_REF][START_REF] Karney | Long-time correlations in the stochastic regime[END_REF][START_REF] Mackay | Stochasticity and Transport in Hamiltonian Systems[END_REF][START_REF] Mackay | Transport in Hamiltonian systems[END_REF][START_REF] Horita | Long-time correlations and expansion-rate spectra of chaos in Hamiltonian Systems[END_REF][START_REF] Ishizaki | Anomalous Diffusion Due to Accelerator Modes in the Standard Map[END_REF][START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF][START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF][START_REF] Stefancich | Anomalous diffusion and ballistic peaks: A quantum perspective[END_REF][START_REF] Kroetz | Hidden high period accelerator modes in a bouncer model[END_REF]). In nonequilibrium statistical mechanics macroscopic properties of matter emerge from microscopic chaotic motion, such as the motion of single atoms or molecules (see e.g. [START_REF] Klages | Microscopic Chaos, and Transport in Nonequilibrium Statistical Mechanics[END_REF] and references therein). Transport analysis in dynamical systems aims to relate the collective motion of ensembles of trajectories between regions of phase space, with several different physical phenomena, such as, the determination of chemical reaction rates, the mixing rates in a fluid, and the particle confinement times in accelerator or fusion plasma devices (see [START_REF]Thirty years of turnstiles and transport, Chaos: An Interdisciplinary[END_REF] for a recent review).

The classical standard map (SM) [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF], also known as the Chirikov map or the kicked rotor system, is a well-known and studied model which exhibits a rather rich diffusive dynamics depending on the value of the so-called kick-strength control parameter K (more information is provided in Sect. 2 below). The system describes the motion of a rotor modeled via the time 1 mgshen002@myuct.ac.za 2 thanos.manos@cyu.fr 3 haris.skokos@uct.ac.za evolution of two variables, respectively its angle and angular momentum, which typically are bounded in the interval [0, 2π) as both variables are considered to be modulo 2π. Depending on the value of the kick-strength K, the angular momentum variable can become 'random' (uncorrelated) performing a Brownian/Gaussian random walk. When chaos dominates the phase space of the system, the motion transport occurs in a fashion termed as 'normal diffusion' in the literature (see e.g. [START_REF] Altmann | Anomalous Transport in Hamiltonian Systems[END_REF]), with the related transport rate referring to the Gaussian distribution of a characteristic quantity of the system, whose variance increases as time evolves.

However, there are intervals of K values where the diffusion deviates from 'normal' and the transport of motion evolves by following an 'anomalous' rate associated with the presence of the so-called accelerator modes (AMs) [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]. The latter ones are regular regions (islands of stability) surrounding stable periodic orbits of different periods in the model's (compact) phase space. There the dynamics produces leaps (equal to 2π or integer multiples of 2π) in the angular momentum variable and orbits get transported infinitely in the de-compactified phase space (i.e. when the modulo requirement is relaxed for the angular momentum). All orbits starting with initial conditions (ICs) inside an AM island are evolving in a ballistic manner, which is linear in time and its diffusion rate is larger than the one observed for normal diffusive transport. Another interme-diate in rate (between normal and ballistic) type of diffusion, which is referred to as superdiffusion, can be observed when ICs happen to lie in the vicinity of the AMs' region or get near them and become trapped as time evolves. In such cases orbits experience the so-called effect of stickiness (see e.g. [START_REF] Dvorak | stickiness" in mappings and dynamical systems[END_REF] and references therein) due to surrounding cantori and exhibit anomalous diffusion. On the other hand, orbits whose ICs lie within ordinary islands of stability do not diffuse and exhibit what is termed as subdiffusion (see e.g. [START_REF] Zaslavsky | Hierarchical structures in the phase space and fractional kinetics: I. Classical systems[END_REF][START_REF] Venegeroles | Leading Pollicott-Ruelle Resonances and Transport in Area-Preserving Maps[END_REF][START_REF] Venegeroles | Leading Pollicott-Ruelle resonances for chaotic areapreserving maps[END_REF][START_REF] Venegeroles | Calculation of Superdiffusion for the Chirikov-Taylor Model[END_REF]).

A rather systematic analysis of the general diffusion processes in the SM for a large set of parameter values generating phase spaces of different diffusion transport processes (normal and anomalous) due to the presence of AMs of different period, was performed in [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF]. That analysis was accompanied by detailed stability maps (where regular and chaotic regions were identified), a description of the momentum distribution with Lévy stable distributions, as well as numerical estimations of the diffusion rates. That work was motivated by studies on the impact of the classical SM on the quantized version of the map and its wave function's localization properties [START_REF] Manos | Dynamical localization in chaotic systems: Spectral statistics and localization measure in the kicked rotator as a paradigm for time-dependent and time-independent systems[END_REF][START_REF] Batistić | The intermediate level statistics in dynamically localized chaotic eigenstates[END_REF], and was extended in [START_REF] Manos | Statistical properties of the localization measure in a finite-dimensional model of the quantum kicked rotator[END_REF] by mainly focusing on quantum (short) characteristic time scales, such as the Heisenberg and localization times.

The global/local transport and the related mean diffusion exponents in the case of a single SM in the presence of AMs were investigated in [START_REF] Harsoula | Global and local diffusion in the standard map[END_REF]. Furthermore, in that work the respective time scales of orbits to escape the initial ballistic motion and converge to normal diffusion, as a function of their distance from the AM island in the phase space, were also studied. That work was later extended in [START_REF] Harsoula | Characteristic times in the standard map[END_REF] where a detailed comparison between the Lyapunov time, the Poincare recurrence time and the escape time was performed. In addition, in [START_REF] Cincotta | Phase correlations in chaotic dynamics: a Shannon entropy measure[END_REF], the authors investigated the presence of phase correlations using the Shannon entropy to distinguish between strong correlations in the time evolution associated with anomalous diffusion, while in [START_REF] Díaz | Controlling escape in the standard map[END_REF] mechanisms for controlling the escape of orbits in the SM were studied. A generalized SM, the so-called rational SM, was considered in [START_REF] Cincotta | Global dynamics and diffusion in the rational standard map[END_REF] (see references therein regarding the model's introduction), where an extensive study and comparison between the different dynamical and diffusion trends appearing in the SM were performed.

Although most works focused on the behavior of a single SM, some studies on higher dimensional maps have also been conducted. For example, processes exhibiting strong chaos and anomalous diffusion in coupled SMs have were studied in [START_REF] Altmann | Hypothesis of strong chaos and anomalous diffusion in coupled symplectic maps[END_REF]. There the existence of an enhanced trapping regime induced by trajectories performing a random walk inside the area corresponding to regular islands of the uncoupled maps was found. Nevertheless, the main focus of that paper was on the dynamics of regions with relatively large size islands of stability, as well as on the determination of diffusion time scales for orbits trapped by invariant tori. In [START_REF] Antonopoulos | Coupled symplectic maps as models for subdiffusive processes in disordered hamiltonian lattices[END_REF] the authors studied the diffusion of weak vs strong chaotic motion in multi-dimensional McMillan coupled maps and demonstrated similarities between such types of symplectic maps and the dynamical properties of disordered Klein-Gordon Hamiltonian systems.

Apart from studies of transport phenomena in systems related to the SM, similar investigations have been performed also for models more directly associated with physical problems. For instance, in [START_REF] Sato | Anomalous diffusion in random dynamical systems[END_REF] a study was carried out on the type of dynamics and anomalous diffusion properties in systems where particles are evolved under the influence of a combined chaotic dynamical system generating Brownian motion-like diffusion and a non-chaotic system in which all particles localize. Recently, complex transport properties have also become relevant in studies of electron motion in low-dimensional nanosystems describing important physical structures, like for example molecular graphene (see e.g. [START_REF] Gil-Gallegos | Energy-dependent diffusion in a soft periodic Lorentz gas[END_REF] and references therein). More specifically, the authors of that work investigated the energydependent diffusion in a soft periodic Lorentz gas with Fermi potentials where macroscopic transport emerges from microscopic chaos and identified different diffusion transport rates at different dynamical regimes.

In this paper, we set out to investigate the long-term (asymptotic) diffusion transport properties of ensembles of ICs in the presence of AMs of different periods in their respective phase space neighborhoods for coupled SMs, extending in some sense the work of [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF]. To this end, we first begin with a rather systematic and exhaustive investigation of the characteristic diffusion trends and properties in the respective uncoupled system (which practically corresponds to a single SM), which possess AMs. We are mainly interested in studying global averaged diffusion rates of the whole phase space and not just diffusion properties of single ICs. We mostly focus on cases with kickstrength values for which the phase space of the single SM is almost globally dominated by chaotic orbits with some small stable islands surrounding AMs. This analysis provides us quite detailed information on the independent dynamics of each map of the coupled system, as well as it allows us to better understand the dissimilarities in the diffusion rates and time-scales that different ensembles of ICs require to converge to their respective long-term (asymptotic) diffusion rate value. Such time scales are found to be a function of the relative fraction of chaotic vs stable areas surrounding an AM.

We then study a system of coupled SMs whose individual SM setups are similar to the ones studied in the first part of our work. Initially, all coupled SMs have the same kick-strength parameter while later on, we also explore SM arrangements of different kick-strengths per map. For these arrangements we investigate (i) the role of the coupling-strength and (ii) the impact of the choice of the ensembles of ICs (different fractions of chaotic regions around the AMs) on each SM, on the global diffusion rates and on the time scales to converge to their longterm values. Finally, we associate these diffusion time scales with the nature of the underlying dynamics, which is quantified through measurements of the extent of the chaotic and regular components. More specifically, we obtain average characterizations of the nature of ICs' ensembles by implementing some well-known and extensively used chaos indicators like the Maximum Lyapunov Exponent (MLE) [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them[END_REF][START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application[END_REF][START_REF] Skokos | The Lyapunov characteristic exponents and their omputation[END_REF] and the Generalized Alignment Index (GALI) method [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method[END_REF][START_REF] Skokos | Detecting chaos, determining the dimensions of tori and predicting slow diffusion in Fermi-Pasta-Ulam lattices by the Generalized Alignment Index Method[END_REF][START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF][START_REF] Skokos | The smaller (SALI) and the generalized (GALI) alignment indices: Efficient methods of chaos detection[END_REF]. These chaos detection techniques have already been successfully used in studies of the chaoticity of coupled SMs [START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF][START_REF] Bountis | Application of the GALI method to localization dynamics in nonlinear systems[END_REF].

The paper is organized as follows: in Sect. 2, we introduce the SM model and the relevant diffusion rate measurements, while in Sect. 3 we provide the definitions of the chaos detection techniques used throughout this work. Our results are presented in Sect. 4 while Sect. 5 summarizes our main findings.

Models and diffusion measures

The two-dimensional (2D) SM is defined by the equations:

x n+1 = x n + y n+1 , y n+1 = y n + K sin(x n ), (1)
where the variables x n and y n are both considered in modulo 2π and n = 0, 1, 2, . . ., counts the number of discrete time steps or in other words the map's iterations. This map represents the evolution of the rotation angle x n and the angular momentum y n of a periodically 'kicked' pendulum rotating in a free field, where the kick by a nonlinear force takes place at each time unit. The constant K measures the kick's nonlinear intensity. The variables x n , y n and the parameter K are dimensionless. Despite the fact that SM is a relatively simple model, analogous to an autonomous Hamiltonian system with 2 degrees of freedom, it exhibits rather rich chaotic and regular dynamics. For K = 0, the SM is linear and as such only quasiperiodic and periodic motion occurs in the system (depending on the initial value y 0 of the angular momentum), filling respectively the 2D phase plane with horizontal straight lines or series of a finite number (equal to the periodicity p of the periodic orbit) of distinct points. For K > 0 some of these lines break into pairs of isolated points, corresponding to periodic orbits, half of which are stable and half unstable [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]. The former are surrounded by families of quasiperiodic orbits (closed curves around these stable periodic orbits), while the latter ones are saddle points surrounded by thin chaotic areas. As the K > 0 value increases, these chaotic areas grow in size and for relatively large K values tend to occupy the whole phase plane.

The quantification of the SM's diffusion process can be done through the relation (see e.g. [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF][START_REF]Thirty years of turnstiles and transport, Chaos: An Interdisciplinary[END_REF]):

(∆y) 2 = D µ (K) n µ , (2)
where ∆y = y ny 0 , with y 0 being the IC of the orbit's angular momentum y n . The quantity (∆y) 2 represents the average of (∆y) 2 at each iteration n, over an ensemble of ICs. The numerically estimated diffusion exponent µ is in the interval [0, 2] and D µ (K), for n → ∞, is the generalized classical diffusion coefficient. When µ = 1 we observe normal diffusion, and D 1 (K) is the so-called normal diffusion coefficient. In cases of anomalous diffusion we encounter two types of behaviors: subdiffusion if 0 < µ < 1 and superdiffusion when 1 < µ ≤ 2, while obviously µ = 0 signifies the lack of diffusion. In the most extreme case where µ = 2, ballistic transport occurs which is strictly associated with the presence and dominance of AMs in the dynamics of the system. The theoretically estimated value of D 1 (K) for normal diffusion (µ = 1) can be found in [START_REF] Izrailev | Simple models of quantum chaos: Spectrum and eigenfunctions[END_REF] to be:

D 1 (K) =        K 2 2 {1 -2J 2 (K) [1 -J 2 (K)]} , if K ≥ 4.5, 0.30(K -K cr ) 3 , if K cr < K < 4.5, (3)
where K cr ≈ 0.9716 and J 2 (K) is a Bessel function. This theoretical result about the dependence of the diffusion coefficient on K fails around the period p = 1 AM existence intervals:

2πl ≤ K ≤ (2πl) 2 + 16, (4)
with l being any positive integer. For example, in the interval obtained from (4) for l = 1 there exist two stable periodic orbits located at y = 0, x = πx * and y = 0, x = π + x * , along with two unstable periodic orbits at y = 0, x = x * and y = 0, x = 2πx * , where x * = arcsin(2π/K).

In order to characterize the different rate of diffusion, we calculate the effective diffusion coefficient [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF][START_REF]Thirty years of turnstiles and transport, Chaos: An Interdisciplinary[END_REF][START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF]]

D eff = (∆y) 2 n , (5)
numerically for large but finite values of n, although the coefficient is theoretically defined for n → ∞. We note that the quantity D eff , in general, is not equal to the coefficient D µ defined in (2) (see also [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF] for more details and references therein). For the 2D SM (1) the dependence of the quantity D eff on the value of K has been already presented in many previous papers, like for example in Figs. 1 of [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF][START_REF] Harsoula | Global and local diffusion in the standard map[END_REF]. In those plots we see peaks in the values of D eff at certain intervals of K values, where AMs are predicted to be present. It is worth noting that these values deviate substantially from those predicted in (3), as the theoretical prediction is based on the assumption of normal diffusion rate.

In our study we also consider a system of N coupled SMs, leading to a 2ND map, described by the following equations [START_REF] Kantz | Internal Arnold diffusion and chaos thresholds in coupled symplectic maps[END_REF]:

x j n+1 = x j n + y j n+1 , y j n+1 = y j n + K j sin(x j n ) -β[sin(x j+1 n -x j n ) + sin(x j-1 n -x j n )], (6)
where j = 1, 2, . . . , N is the index of each SM, K j its respective nonlinear kick-strength parameter, and β the coupling-strength parameter between neighboring maps. The choice of each K j value plays a pivotal role since it initializes each map at the dynamical state of the respective 2D SM (1) before the coupling takes place. This in turn is related to the initial amount of regular and chaotic motion, as well as the presence of AMs in its respective phase plane. In our investigations, we consider periodic boundary conditions, i.e. x 0 n = x N n and x N+1 n = x 1 n , a choice which preserves the symplectic nature of the coupled system.

In order to investigate the diffusion process for the N coupled SMs, we use a generalized expression similar to Eq. ( 5), defining (as a limit for n → ∞) the effective diffusion coefficient of the 2ND map as:

D N eff = 1 n N j=1 (y j n+1 -y j 0 ) 2 = 1 n N j=1 (∆y j ) 2 . ( 7
)
In this expression the sum is over the total number of maps and each map has the same number of ICs, but the location of each ensemble of ICs may vary between the 2D SMs.

Chaos measure techniques

Let us now discuss the main numerical techniques we implement in our study, starting with a brief presentation of the two chaos detection methods we use, namely the MLE and the GALI, emphasizing their application to area-preserving, symplectic maps. In our presentation we refer to the more general case of the 2ND map [START_REF] Karney | Long-time correlations in the stochastic regime[END_REF], which, for simplicity, can be expressed as:

x n+1 = F(x n ), (8)
with x n denoting the vector of the map's coordinates at iteration n and F(x n ) the functional form of the map, i.e. the expressions at the right hand side of [START_REF] Karney | Long-time correlations in the stochastic regime[END_REF].

The regular or chaotic nature of the map's orbits is determined by the evolution of small perturbations to these orbits. Such perturbation defines the so-called deviation vector w n . The propagation of this vector is governed by the system's tangent map:

w n+1 = ∂F ∂x (x n ) • w n . (9)
The MLE λ is the average expansion (or shrinking) rate of the distance between two neighboring orbits in the phase space of a dynamical system (quantified in a first order approximation by the length of the deviation vector) and it is computed as [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them[END_REF][START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application[END_REF][START_REF] Skokos | The Lyapunov characteristic exponents and their omputation[END_REF]:

λ = lim n→∞ Λ(n), (10)
where Λ(n) is the so-called finite-time MLE (ftMLE) given by:

Λ(n) = 1 n ln w n w 0 , (11)
with w 0 and w n being respectively deviation vectors from the given orbit at iterations n = 0 and n > 0, and • denoting any vector norm (in our study we implement the usual Euclidean norm). In the case of chaotic orbits λ is positive, while it is zero for regular orbits. Since it is impractical to actually compute λ (10), as it is obtained through a limit when the number of iterations n goes to infinity, we rely on computations of the ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF]. In practice, for chaotic orbits Λ(n) tends to a positive value, while it goes to zero as Λ(n) ∝ ln n n ∼ n -1 for regular orbits.

Although the MLE has been extensively used as a chaos indicator one practical problem it often faces is the potential slow convergence of Λ(n) to its limiting λ value. Over the years several methods which overcome this problem have been developed like the Smaller Alignment Index (SALI) [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF][START_REF] Skokos | How does the Smaller Alignment Index (SALI) distinguish order from chaos?[END_REF][START_REF] Skokos | Detecting order and chaos in Hamiltonian systems by the SALI method[END_REF] and its generalization the GALI, the Fast Lyapunov Indicator (FLI) [START_REF] Froeschlé | Fast Lyapunov indicators. Application to asteroidal motion[END_REF][START_REF] Froeschlé | The fast Lyapunov indicator: a simple tool to detect weak chaos. Application to the structure of the main asteroidal belt[END_REF][START_REF] Froeschlé | On the structure of symplectic mappings. The Fast Lyapunov Indicator: a Very Sensitive Tool[END_REF][START_REF] Lega | Theory and applications of the fast Lyapunov indicator (FLI) method[END_REF] and its variants [START_REF] Barrio | Sensitivity tools vs. poincaré sections[END_REF][START_REF] Barrio | Painting chaos: a gallery of sensitivity plots of classical problems[END_REF][START_REF] Barrio | Theory and applications of the orthogonal fast lyapunov indicator (ofli and ofli2) methods[END_REF], the mean exponential growth of nearby orbits (MEGNO) [START_REF] Cincotta | Simple tools to study global dynamics in nonaxisymmetric galactic potentials -I[END_REF][START_REF] Cincotta | Phase space structure of multidimensional systems by means of the mean exponential growth factor of nearby orbits[END_REF][START_REF] Cincotta | Theory and applications of the mean exponential growth factor of nearby orbits (MEGNO) method[END_REF] and the '0-1' test [START_REF] Gottwald | A new test for chaos in deterministic systems[END_REF][START_REF] Gottwald | Testing for chaos in deterministic systems with noise[END_REF][START_REF] Gottwald | The 0-1 test for chaos: A re-view[END_REF] to name a few (for a collection of review papers on various modern chaos detection techniques the reader is referred to [START_REF]Chaos Detection and Predictability[END_REF]). Here we will implement the GALI method, which has been proven to be very efficient and has already been successfully applied in studies of various dynamical systems (see e.g. [START_REF] Bountis | Application of the GALI method to localization dynamics in nonlinear systems[END_REF][START_REF] Manos | Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models[END_REF][START_REF] Manos | Scaling with system size of the Lyapunov exponents for the Hamiltonian mean field model[END_REF][START_REF] Manos | Interplay between chaotic and regular motion in a time-dependent barred galaxy model[END_REF][START_REF] Chaves-Velasquez | Boxy orbital structures in rotating bar models[END_REF]).

For 2ND maps the GALI of order k (GALI k ), 2 ≤ k ≤ 2N, is determined through the evolution of k initially linearly independent deviation vectors w i 0 , i = 1, 2, . . . , k. To avoid overflow problems, the resulting deviation vectors w i n , i = 1, 2, . . . , k, are continuously normalized, but their directions are kept intact. Then, according to [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method[END_REF] GALI k is defined as the volume of the k-parallelogram having as edges the k unit deviation vectors ŵi n = w i n / w i n , i = 1, 2, . . . , k, determined through the wedge product of these vectors as:

GALI k (n) = ŵ1 n ∧ ŵ2 n ∧ • • • ∧ ŵk n , (12)
with • denoting the usual Euclidean norm. From this definition it is evident that if at least two of the deviation vectors become linearly dependent, the wedge product in [START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF] becomes zero and GALI k vanishes. As mentioned earlier, the GALI method is a generalization of the SALI technique, which practically coincides with GALI 2 [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method[END_REF]. The behavior of the GALI k for regular and chaotic orbits was theoretically studied in [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method[END_REF][START_REF] Skokos | Detecting chaos, determining the dimensions of tori and predicting slow diffusion in Fermi-Pasta-Ulam lattices by the Generalized Alignment Index Method[END_REF], where it was shown that all GALI k (n) tend exponentially to zero for chaotic orbits, with exponents which depend on the first k Lyapunov exponents of the orbit. On the other hand, in the case of regular orbits, GALI k (n) remains practically constant and positive if k is smaller or equal to the dimensionality of the torus on which the motion occurs, otherwise, it decreases to zero following a well-defined power law decay.

In the particular case of 2D maps, like the SM (1), GALI 2 (which is the only GALI that can be defined in this case) tends to zero both for regular and for chaotic orbits, following however completely different time rates, which allow us to distinguish between the two cases [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF]. More specifically, in the case of a chaotic orbit, any two deviation vectors will be aligned to the direction defined by the MLE λ, and consequently GALI 2 will eventually tend to zero following an exponential decay of the form GALI 2 (n) ∝ e -2λn . In the case of regular orbits, any two deviation vectors tend to fall on the tangent space of the torus on which the motion lies. For a 2D map, this torus is a 1D invariant curve, whose tangent space is also 1D and consequently any two deviation vectors will become aligned. Thus, even in the case of regular orbits in 2D maps GALI 2 tends to zero. This decay follows the power law GALI 2 (n) ∝ 1/n 2 [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF][START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF]. It is worth noting that, similarly to chaotic orbits, the GALI 2 (n) tends exponentially to zero for unstable periodic orbits, while it remains practically constant, fluctuating around a positive value, for stable periodic orbits [START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF].

The behavior of the GALI 2 for all these different types of orbits is seen in Figs. 1(b) and (e). More specifically in Fig. 1(a) the phase space portrait of SM (1) for K = 3.1 is depicted, where a chaotic (red points) and a regular (blue points) orbit are also identified. The evolution of the GALI 2 for these orbits is shown in Fig. 1(b) with the index exhibiting an exponential decay for the chaotic orbit (red curve) and a power law decrease proportional to n -2 (indicated by the dashed line) for the regular orbit (blue curve). From the results of Fig. 1(c) we see that the ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] tends to zero for the regular orbit (blue curve), eventually following a Λ(n) ∝ n -1 law (indicated by the dotted line), while it saturates to a positive value for the chaotic d), (e) and (f)] K = 6.5. In (a) and (d) phase space portraits of the SM are shown, which are produced by 10 6 iterations of several ICs (black points), while colored points are used to identify particular orbits. In (a) the consequents of the chaotic orbit with IC (x 0 , y 0 ) = (2, 1.5) (red pentagon point) are plotted in red, while blue points correspond to a regular orbit with IC (x 0 , y 0 ) = (3, 0.5) (blue diamond point). In (d) the location of two AMs of period p = 1 are identified, a stable AM (green circle point) at (x 0 , y 0 ) = (1.8298, 0) and an unstable AM (purple square point) at (x 0 , y 0 ) = (1.3118, 0). The evolution of the GALI 2 of all these orbits, with respect to the number of iterations n, is illustrated in (b) for the chaotic (red curve) and the regular (blue curve) orbit of (a), and in (e) for the stable (green curve) and the unstable AM (purple curve) of (d). In (b) GALI 2 (n) tends to zero exponentially fast for the chaotic orbit, while it decreases to zero following a power law GALI 2 ∝ n -2 (indicated by a dashed line) for the regular orbit. In (e) GALI 2 (n) remains practically constant for the stable AM, while it tends to zero exponentially fast for the unstable AM. In (c) and (f) we see the evolution of the ftMLE Λ(n) [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF], respectively for the chaotic (red curve) and the regular (blue curve) orbit of (a), and for the stable (green curve) and the unstable (purple curve) AM of (d). Λ(n) saturates to a positive value in (c) and (f), respectively for the chaotic orbit and the unstable AM, while it tends to zero following a power law Λ(n) ∝ n -1 (indicated by a dotted line in both panels) for the regular orbit and the stable AM. orbit (red curve). From the comparison of Figs. 1(b) and(c) it is worth noting how much faster the evolution of the GALI 2 allows the clear distinction between the regular and chaotic behavior of orbits, with respect to the ftMLE. More specifically, we see that after n ≈ 50 iterations the GALI 2 of the chaotic orbit becomes practically zero (GALI 2 ≈ 10 -15 ) attaining values several orders of magnitude smaller than in the case of the regular orbit, while at the same time a clear distinction between the two orbits cannot be securely established from the results of Fig. 1(c), as for example the ftMLE of the regular orbit has not manifested a clear and well-defined tendency to decrease.

In Fig. 1(d) the phase space location of a stable (green circle point) and an unstable (purple square point) AM for the 2D SM (1) with K = 6.5 is shown. In accordance to [START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF] the GALI 2 of the stable AM remains practically constant [green curve in Fig. 1(e)], while for the unstable AM [purple curve in Fig. 1(e)] it decreases exponentially fast to zero, showing a similar behavior to the one of the chaotic orbit in Fig. 1(b), as in the vicinity of the unstable AM chaotic motion occurs. In Fig. 1(f) we see that the ftMLE Λ for the stable (green curve) and the unstable (purple curve) AM behaves respectively as in the case of the regular and chaotic orbit of Fig. 1(c). We note that for the calculation of the GALI 2 and the ftMLE we impose the modulo 2π requirement in order to avoid numerical problems due to the divergence of orbits to infinity.

Results

In this study, we set out to ultimately investigate the rate of diffusion and the respective dynamics of ensembles of orbits in coupled SMs. However, and before doing that, we will perform an extensive study of the dynamics of the single SM (1) focusing on the effects of AMs (of different periods p) in the long term dynamics and transport properties of various sets of ICs.

Single standard map

In order to perform the numerical evaluation of the diffusion exponent µ of (2) for a specific region of the SM's phase space we consider a grid of 315 × 315 equally spaced ICs (i.e. ≈ 100, 000 ICs) and calculate the evolution of the variance (∆y) 2 [Eq. ( 2)] as a function of the map's iterations n. In Fig. 2(a) we present the outcome of this approach for sets of orbits with ICs in small regions around the four particular orbits considered in Fig. 1. More specifically we perform this computation for a small neighborhood of the regular (blue curve) and the chaotic orbit (red curve) of Fig. 1(a), as well as for regions containing the IC of the unstable (purple curve) and the stable AM (green curve) of Fig. 1(d).

From the results of Fig. 2(a) we see that ballistic transport 2) of the angular momentum y n of the 2D SM (1) as a function of the map's iterations n, for orbits in four different small regions of the map's phase space. For each region a set of ICs on a 315 × 315 grid of equally spaced mesh points (i.e. in total ≈ 100, 000 ICs) are taken and the variance is obtained as average over all these orbits. For K = 3.1 ICs in the region appears in the region around the stable AM as the increase of (∆y) 2 (green curve) is well described by a function proportional to n 2 (dashed line), i.e. having a diffusion exponent µ = 2. On the other hand, the dynamics around the unstable AM (purple curve) and the chaotic orbit (red curve) practically leads to normal diffusion (µ = 1) as both curves seem to follow quite well a (∆y) 2 ∝ n increase represented by a dotted line, while no diffusion is observed in the immediate neighborhood of the considered regular orbit (blue curve). Note that typical regular orbits and stable AMs have similar properties from a dynamical point of view, see e.g., the evolution of the ftMLEs given by blue and purple curves in Figs. 1(c) and (f), however their diffusion rates are different as are respectively presented in Figs. 2(a) and (b). The former ones do not practically diffuse (µ ≈ 0) while the latter ones diffuse with extreme rate values (µ ≈ 2).

The value of µ can be numerically estimated by performing a linear fitting of the logarithm of (∆y) 2 with respect to log 10 n, but in order to identify potential changes in the dynamical evolution of the ensembles of ICs, which will be reflected in the value of the diffusion exponent, we prefer to monitor the evolution of the µ value as n increases by first smoothing the variance of the momentum values through a locally weighted regression smoothing algorithm [START_REF] Cleveland | Locally weighted regression: An approach to regression analysis by local fitting[END_REF]. Then by implementing the numerical approach followed in [START_REF] Laptyeva | The crossover from strong to weak chaos for nonlinear waves in disordered systems[END_REF][START_REF] Bodyfelt | Nonlinear waves in disordered chains: Probing the limits of chaos and spreading[END_REF] we compute µ in a running window in log 10 n, whose width is much smaller than the total length of the system's evolution. Implementing the latter approach to the results presented in Fig. 2(a) we create Fig. 2(b), where we clearly observe the presence of ballistic transport (µ = 2, green curve) and practically normal diffusion (µ ≈ 1, purple and red curves), as well as the absence of it (µ = 0, blue curve). From the presented results we see that in all cases no significant changes in the evolution of µ are observed, because after some transient phase up to n ≈ 10 3 (apart from the ballistic case where µ remains constant from the beginning) the values of µ do not notably vary, as can also be deduced from the results of Fig. 2(a). Nevertheless, as we will see below [e.g. in Fig. 4(b) and the last column of Fig. 7] plots similar to Fig. 2(b) will be extremely helpful in identifying the possible alterations in the diffusion and transport properties of sets of ICs, as well as in determining the number of iterations needed for µ to approach its asymptotic µ = 2 value.

Extending the analysis presented in [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF] we compute the dependence of the diffusion exponent µ on the value of K for a set of 315 × 315 equally spaced ICs on the entire phase space of the SM, i.e. x, y ∈ [0, 2π), for many more iterations than the ones considered in [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF] (actually up to n = 10 7 ), in order to also investigate the dependence of the value of µ on n. In Figs. 3(a) and (c) we present µ as a function of K when, respectively, n = 10 4 and n = 10 5 iterations of the considered set of orbits are performed. In particular, µ is estimated through a linear fitting of log 10 (∆y) 2 with respect to log 10 n in the last two decades of the evolution. This means that, for example in the case of Fig. 3(c) the fitting is performed for data with 3 ≤ log 10 (n) ≤ 5.

At the K intervals [Eq. ( 4)] where AMs of period p = 1 are expected [also indicated at the horizontal black line at the bottom of Figs. 3(a) and(c)] we observe high peaks denoting anomalous diffusion due to the presence of AMs. We note that for the considered values of K ≤ 70, 11 such peaks appear. The largest obtained values µ * at these peaks are marked by black filled circles. The direct comparison of Figs. 3(a) and(c) shows that these peaks, and the related µ * values, increase as the number of performed iterations grows, approaching the µ = 2 value, which corresponds to ballistic transport. We note that at the same intervals of K values the computed diffusion coefficient D eff (5) also exhibits abrupt spikes [Figs. 3(b) and (d)], whose height increases with n as well. The other peaks observed in Figs. 3(a) and (c) for K < 4π (whose height also increases for larger n) are related to AM of periods p > 1, while K intervals with µ ≈ 1 are associated with normal diffusion processes. We emphasize here that the ballistic motion caused by AM ICs prevails on the system's global diffusion and the effective diffu- sion coefficient D eff (5) will go to infinity, as was also reported in [START_REF] Harsoula | Global and local diffusion in the standard map[END_REF], while the corresponding, entire phase space diffusion exponent µ (2) will eventually converge to the value 2. This is valid even for very large K values, as the AM islands exist even when K → ∞, having nevertheless smaller and smaller sizes as K grows.

In order to further investigate the dynamical properties of the SM when anomalous diffusion is observed due to the presence of period p = 1 AMs, we study in more detail the behavior of the map for the 11 K values at which the peaks in Fig. 3 appear. In Fig. 4 we present results for only 4 such cases [namely for K = 6.866 (blue curves), K = 25.291 (green curves), K = 44.002 (red curves) and K = 69.173 (black curves)] in order to not overload the plots. More specifically, the dependence of (∆y) 2 and the related diffusion exponent µ [Eq. ( 2)] on the number n of map's iterations is respectively shown in Figs. 4(a) and (b). From these results we see that in all cases µ eventually becomes µ = 2 (indicating asymptotic ballistic transport), but the number of iterations needed to reach this value increases as K becomes larger. This feature is clearly seen in Fig. 4(b) as the presented curves saturate to µ = 2 at larger n for higher K values.

It is also worth noting that the initial phases of the dynamical evolution are better described as normal diffusion when K has large values (green and in particular red and black curves in both panels of Fig. 4) as for these cases the increase of (∆y) 2 is practically proportional to n [dotted line in Fig. 4(a)] and the corresponding µ values are close to µ = 1 [Fig. 4(b)]. This delayed kick-in of anomalous diffusion for larger K values is also reflected on the lowering of the peak heights in Figs. 3(a) and (c) when K increases.

The increase of the peak highest values µ * of the diffusion exponent in the first 11 intervals [Eq. ( 4)] where p = 1 AMs exist, and their tendency to approach µ = 2 as n grows, is clearly seen in Fig. 5(a) where results for n = 10 3 (blue diamonds), n = 10 4 (green pentagons), n = 10 5 (red squares), n = 10 6 (orange triangles) and n = 10 7 (purple circles) are plotted. We note that the points for n = 10 4 and n = 10 5 are the ones respectively marked by black circles in Figs. 3(a) and(c). The solid curves correspond to fittings of each set of data points by a power law of the form:

µ * = AK B , (13)
and the dependence of the fitting parameters A and B on the iterations n of the map is respectively seen in Figs. 5(b) and (c).

From the results of Fig. 5(a) we vividly see that the highest µ values in each peak tend asymptotically to µ = 2 as n increases. This feature is also reflected in the values of parameters A and B, which respectively approach A = 2 [dashed line in Fig. 5(b)] and B = 0 [dashed line in Fig. 5(c)]. These tendencies denote that for large n the values of µ * actually become independent of K and almost equal to µ = 2 as (13) gives µ * = 2 for A = 2 and B = 0. Nevertheless, as indicated by the results of Figs. 4(b) and 5(a) the diffusion exponents associated with the larger K peaks require more iterations to reach the asymptotic, limiting value µ = 2. This behavior is clearly seen in Fig. 6(a) where we plot the number n of iterations needed for µ to reach µ = 2 (actually n is defined as the number of iterations needed for µ to take values very close to µ = 2 and in particular to become µ ≥ 1.992). In the inset of Fig. 6(a) we plot n as a function of the length ∆K l = (2πl) 2 + 16 -2πl of the period p = 1 AM existence intervals (4), and since ∆K l decreases for larger K (larger l), n values are higher for smaller ∆K l .

In order to relate the delay of the diffusion exponent to reach the value µ = 2 with the strength of chaos of the SM we compute the average Λ ftMLE (11) over a set of orbits covering the entire phase space, as an overall indicator of the system's chaoticity, for the 11 K values considered in Fig. 5 space, i.e. x, y ∈ [0, 2π), for n = 10 5 iterations, registering the average value Λ of the ftMLE at the end of the simulations, considering also one standard deviation of this process as a good indicator of the error of the computed average. The output of this procedure is seen in Fig. 6(b) where a clear increase of Λ as K grows is evident. Actually, the obtained results are in excellent agreement with the relation Λ = ln(K/2) [dashed red curve in Fig. 6(b)] which has been presented in previous works [START_REF] Shevchenko | On the maximum Lyapunov exponent of the motion in a chaotic layer[END_REF][START_REF] Harsoula | Characteristic times in the standard map[END_REF]. Thus, we conclude that in general the existence of period p = 1 AMs eventually leads to ballistic transport (µ = 2) for sets of orbits distributed in the entire phase space of the SM but for stronger chaoticities (larger K values) this behavior is achieved after longer iteration intervals. We note that in the computation of the average Λ value we potentially mix regular and chaotic orbits. Nevertheless, although islands of stability always exist in the phase space of the SM, their extent decreases rapidly as K grows and the overall dynamics is practically defined by chaotic motion (see e.g. [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]). Thus, the possible inclusion of some regular orbits, whose ftMLE tends to zero, does not practically influence the computation of Λ , as the used ensemble of orbits (10, 000) is sufficiently large. In addition, the number of iteration used (n = 10 5 ) for the evaluation of Λ is large enough to permit a reliable estimation of a representative value for the system's chaoticity, as Λ has already saturated to its limiting value for all considered cases [see inset of Fig. 6(b)].

Let us now study in more detail the relation between the SM's chaoticity and its diffusion properties in the presence of AMs of various periods. The different behavior of the GALI 2 for regular and chaotic orbits, seen in Fig. 1(b), allows us to use this index for a fast and clear distinction between regions of chaos and order in the 2D phase space of the SM (1). More specifically, in order to estimate the percentage of chaotic orbits of the SM for a given value of K we implement the subsequent approach: We follow the evolution of orbits whose ICs lie on a 2D grid with equally spaced points on the 2D phase space of the map, and register for each orbit the value of GALI 2 after n = 50 iterations. All orbits having values of GALI 2 smaller than 10 -8 at n = 50 are characterized as chaotic, while all others are considered as non-chaotic (we note that a similar process was also performed in [START_REF] Manos | Survey on the role of accelerator modes for anomalous diffusion: The case of the standard map[END_REF]). We also point out that for each orbit a different random set of orthonormal initial deviation vectors is used setting GALI 2 (0) = 1. The 10 -8 threshold value has been broadly and successfully used in previous works (see for example [START_REF] Manos | Studying the global dynamics of conservative dynamical systems using the SALI chaos detection method[END_REF][START_REF] Manos | Application of the generalized alignment index (GALI) method to the dynamics of multi-dimensional symplectic maps[END_REF][START_REF] Manos | Global dynamics of coupled standard maps[END_REF]), and can safely distinguish the two types of orbits as is also evident from the results of Fig. 1(b). In addition, we have checked that n = 50 iterations are sufficient to securely distinguish between the exponential (chaotic motion) and the power law (regular motion) decay of GALI 2 .

For the single SM (1) we start by choosing a K value for which AMs are present and use the GALI 2 method to detect the regular and chaotic regions in the system's phase space, producing color maps which clearly give a global overview of the dynamics and measuring at the same time the fraction P C of regions where chaotic motion occurs. In order to create such detailed color plots we consider ensembles consisting of 250, 000 ICs on a grid of 500 × 500 cells (i.e. one orbit per cell) covering the entire phase space, i.e. (x, y) ∈ [0, 2π) × [0, 2π), or a subspace of it. The evolution of the orbits of this grid is also used to estimate the percentage P C of chaotic orbits through the computation of their GALIs. The same dense grid is also used to measure the 'local' rate of diffusion for each cell and create color plots of the related µ value [Eq. ( 2)], but for such computations we additionally consider 50 × 50 ICs inside each one of the 500 × 500 cells (i.e. 2,500 orbits per cell) and numerically calculate the diffusion exponent µ for each cell separately by using (2) and performing a linear fit of the variance (∆y) 2 as a function of the iterations n in log-log scale for 3 ≤ log 10 (n) ≤ 4. These processes produce detailed pictures, with fine resolutions, allowing the systematic analysis of both measurements (the GALI index for the system's chaoticity and the diffusion exponent µ for the diffusion properties of the map) within relatively feasible CPU times and they will be followed whenever we create such color plots in this work. In this way we are able to produce comprehensive plots where regular and chaotic regions are clearly indicated, as well as detect regions characterized by different diffusion rates.

Apart from the global study of the 1D SM we also focus our attention on subspaces of the entire phase space, and more specifically on regions which contain AMs (where superdiffusion takes place) surrounded by chaotic areas (where normal diffusion occurs). By gradually changing the size of the considered subspace and consequently the percentage P C of chaos (i.e. the extent of the chaotic region) around the AMs we investigate the role and impact of this ratio on the diffusion transport properties of the considered subspaces.

The outcomes of these computations for four K values corresponding to the presence of AMs of periods p = 1, 2, 3 and 4 (respectively for K = 6.5, K = 4.0844, K = 6.9115 and K = 3.1) are shown in Fig. 7. In each case the first two panels correspond to color plots of phase space regions based on the value of GALI 2 , with the next two panels depicting the same regions colored according to the computed µ values. In the GALI 2 color plots, regions of regular behavior correspond to large GALI 2 values colored in yellow, while dark purple and black areas denote chaotic behavior. In the µ color plots, regions where ballistic diffusion takes place (µ ≈ 2) are colored in yellow, while dark red and purple colored areas respectively correspond to normal diffusion (µ ≈ 1) and subdiffusing processes.

For each K case the position of the stable AM is denoted by a small black rectangle point in the right GALI 2 and µ color plots (respectively second and fourth column in Fig. 7). Examples of AMs of different periods, as well as the presentation of a method to locate them can be found in [START_REF] Harsoula | Global and local diffusion in the standard map[END_REF].

For each K case we consider 6 phase space regions whose positions are shown in the GALI 2 and µ color plots of Fig. 7 by different rectangles. The x and y ranges of these rectangles, along with their names are listed in Table 1. The names of these regions are composed by the value p of the AM's period followed by a letter from 'A' to 'F'. The various regions were chosen so that the percentage P C of chaotic orbits they contain decreases as we move from 'A' to 'F', taking care at the same time that regions with similar names have equivalent P C values, namely P C ≈ 99% (A), 75% (B), 50% (C), 25% (D), 12.5% (E), 6.5% (F) (the exact P C values are also reported in Table 1). For example, for K = 6.5 (upper row in Fig. 7) regions 1A, 1B and 1C are seen in the left GALI 2 and µ color plots, where they are respectively denoted by dark pink, light blue and black bordered rectangles. Regions 1D, 1E and 1F are respectively depicted in the right GALI 2 and µ panels by green, light red and cyan rectangles.

The panels in the rightmost column of Fig. 7 show the diffusion exponent µ [Eq. ( 2)] as a function of the number of iterations n for every considered phase space region of each K case. In all these panels results for similar P C values are presented by curves with the same color. More specifically, dark pink, light blue, black, green, light red and cyan curves respectively correspond to 'A', 'B', 'C', 'D', 'E' and 'F' named regions. From these results it becomes again clear that the diffusion exponent µ of ensembles containing more chaotic orbits (larger P C values) requires more iterations n to approach its limiting value µ = 2. The corresponding n values are reported in the last column of Table 1.

We note here that the GALI 2 color plots in Fig. 7 were created for n = 50 iterations, which, as we have already argued, are sufficient to distinguish between regular and chaotic orbits, while the µ color plots are obtained by the analysis of data for 10 3 ≤ n ≤ 10 4 . Thus, both types of plots capture the system's dynamical features for these specific number of iterations. As n increases the dynamics evolves and some structures in these color plots will alter. For example, the red colored points inside the purple colored (chaotic) regions of the GALI 2 plots, as well as the purple points inside the red colored areas of the µ plots correspond to stickiness effects along unstable asymptotic curves emanating from unstable periodic orbits around the island of stability (the so-called stickiness in chaos phenomenon [START_REF] Contopoulos | Stickiness in chaos[END_REF][START_REF] Contopoulos | Stickiness effects in chaos[END_REF]), but for a very large number of iterations they will clearly reveal their chaotic nature obtaining very small GALI 2 Table 1: The number n of iterations of the 2D SM (1) needed for the diffusion exponent µ [Eq. ( 2)] to practically reach its maximum asymptotic µ = 2 value (ballistic transport) for several K values and ensembles of ICs in different phase space regions containing various percentages P C of chaotic orbits. The considered K values are K = 6.5, K = 4.0844, K = 6.9115, K = 3.1 and respectively correspond to the appearance of stable AMs of period p = 1, 2, 3 and 4. The names of the various regions are composed by the value p of the AM's period followed by a letter from 'A' to 'F'. The ranges in the x and y coordinates for each region were chosen so that ensembles with the same letter have practically the same P C value, which decreases as we move from 'A' to 'F' [P C ≈ 99% (A), ≈ 75% (B), ≈ 50% (C), ≈ 25% (D), ≈ 12.5% (E), ≈ 6.5% (F)].

K

Ensemble name: 12) value (panels in the two left columns) after n = 50 iterations of orbits with ICs on a dense grid covering the considered phase space region, and the diffusion exponent µ [Eq. ( 2)] (panels in the two middle columns) computed by appropriately fitting the variance (∆y) 2 (2) values for iterations n in the interval 10 3 ≤ n ≤ 10 4 (see text for more details). Each row of panels corresponds to a different K value, namely K = 6.5, K = 4.0844, K = 6.9115 and K = 3.1 (from top to bottom), for which stable AMs of period p = 1, 2, 3 and 4 respectively exist. Points are colored according to the color scales on top of the panels. In the GALI 2 plots yellow regions refer to regular motion, dark blue and black areas correspond to chaotic motion, while in-between colors to weakly chaotic ICs requiring more iteration to reveal their chaotic nature. In the µ color plots ICs leading to subdiffusion rates are colored in dark blue (0 ≤ µ < 1), normal diffusion cases are depicted in dark red (µ ≈ 1), superdiffusion corresponds to lighter shades of red (1 ≤ µ < 2) and ballistic transport is denoted by yellow (µ ≈ 2). The positions of the stable p = 1 [(1.8298, 0)], p = 2 [(4.211, 4.9324)], p = 3 [(1.80, 0.0)] and p = 4 [(0.3486, 2.6121)] AMs are denoted by small black rectangular points in the right GALI 2 and µ panels. In each case 6 different regions of the phase space are considered. The names of these regions, along with their ranges in the x and y coordinates, are given in Table 1, and their positions are indicated by different rectangles in the color plots accompanied by the names of the regions. These regions are ordered in Table 1 in decreasing content of chaotic orbits (quantified by the percentage P C of ICs leading to chaotic motion) and they were chosen so that in all cases ensembles with similar names have equivalent P C values, i.e. P C ≈ 99% for 'A' regions, P C ≈ 75% for 'B', P C ≈ 50% for 'C', P C ≈ 25% for 'D', P C ≈ 12.5% for 'E' and P C ≈ 6.5% for 'F'. The panels in the last column depict the diffusion exponent µ [Eq. ( 2)] of the 6 different regions as a function of the map's iterations n.

x
values, also resulting to µ ≈ 1. More generally, for large n GALI 2 plots will practically be covered by two colors: yellow (regular motion) and black (chaotic motion), while in the µ plots we will actually be able to see only regions leading to ballistic (yellow points, µ ≈ 2) or normal (red points, µ ≈ 1) diffusion.

Coupled standard maps

So far, we have performed a detailed and systematic investigation of diffusion trends, characteristic time scales and properties of single 2D SMs, which possess AMs of different periods, using ensembles of ICs which are dominated by the presence of chaotic orbits (of various fractions) around these AMs. Now we set out to investigate the way that the average diffusion rate of sets of ICs and the average chaoticity of a system of coupled SMs ( 6) is affected by the presence of AMs. In more detail, we aim to better understand how the diffusion and chaoticity are influenced by the individual dynamics of each coupled map (defined by the values of the respective kick-strength parameters, K j , j = 1, 2, . . . , N), coupling-strength (quantified by the value of the β parameter), as well as different maps' arrangements. To this end, we consider a relatively small in size (N = 5) system of coupled 2D SMs with different setups, i.e. uniformly equal kick-strengths for all maps, or maps with different kickstrengths. This choice of the system's size allows us to study sufficiently large sets of ICs, which are needed for the reliable diffusion rate estimation, within feasible CPU times.

For the system of N = 5 coupled SMs (6) we, in general, consider ensembles of ICs whose projections in each one of the N 2D SMs are on a 315×315 equally spaced grid (i.e. ≈100,000 for each one of the coupled 2D SMs). Due to the type of the interaction between neighboring maps [see Eq. ( 6)] we pay special attention in avoiding the zeroing of the coupling interaction, which appears by default whenever the angle coordinates x j 0 , j = 1, 2, . . . , N, of the IC are set to be equal in neighboring 2D SMs, i.e. whenever x j 0 = x j+1 0 . Several strategies can be implemented to avoid this situation whenever the same range [x min , x max ] of x j 0 values is considered in adjacent maps. In our study, we impose a shift in the angular position variable x of the form x j+1 0 = x j 0 + d, where d = (x maxx min )/N is a displacement defined by dividing the angular value range of the chosen ensemble of ICs by the total number N of maps in the system. If the computed x j+1 0 value falls outside the interval [x min , x max ], i.e. x j+1 0 > x max , then this value is reduced by (x maxx min ) so that it is moved inside the appropriate value range.

For the numerical estimation of the generalized diffusion exponent µ we use a similar expression to (2) (which is valid for a single 2D SM), where the left-hand side of the equation is replaced by the quantity N j=1 (∆y j ) 2 appearing in [START_REF] Mackay | Stochasticity and Transport in Hamiltonian Systems[END_REF], i.e.

N j=1

(∆y j ) 2 = D N µ n µ , (14)
with D N µ being (for n → ∞) a diffusion coefficient analogous to D µ in [START_REF] Rechester | Calculation of turbulent diffusion for the Chirikov-Taylor model[END_REF]. We also remark that we have checked that the considered sizes of ensembles of ICs used in our investigations manage to correctly capture the systems' dynamics in feasible CPU times.

Coupling identical 2D standard maps

Following a similar approach as for single SMs (see Fig. 3), we start by calculating the diffusion exponent µ [Eq. ( 14)] and the effective diffusion coefficient D N eff [START_REF] Mackay | Stochasticity and Transport in Hamiltonian Systems[END_REF] for the coupled system of SMs. Our aim is to detect conditions and/or parameter values under which one may expect to observe global diffusion properties similar to those exhibited for single SMs. In other words, we want to investigate the impact of the coupling between SMs in the long-term diffusion properties of the system's phase space, as well as the time (iteration) intervals for which the coupled system may still be influenced by the presence of AMs in one or more of the single 2D SMs.

We first consider a system of coupled SMs with similar individual setups for each one of its 2D maps, i.e. keeping equal kick-strength parameters for all maps to K j = K, j = 1, 2, . . . , 5. In Fig. 8(a) and in its inset, we respectively plot the diffusion exponent µ [Eq. ( 14)] and the effective diffusion coefficient D N eff [START_REF] Mackay | Stochasticity and Transport in Hamiltonian Systems[END_REF] as functions of K, for map [START_REF] Karney | Long-time correlations in the stochastic regime[END_REF] with N = 5 and β = 10 -4 (a moderate coupling value allowing dynamical effects to take place within feasibly CPU times). These results are obtained for a set of ICs on a 315 × 315 grid covering the whole phase space of each one of the coupled 2D SMs. We note that, after experimenting with different ensemble sizes (total number of ICs), we found that the considered number of ICs is sufficiently large to correctly capture the basic dynamical features of the system, and at the same time small enough to allow the performance of extensive numerical computations in feasible computational times. Both the diffusion exponent µ and the effective diffusion coefficient D N eff are computed after n = 10 4 iterations. Fig. 8(b) depicts a similar analysis to the one of Fig. 8(a) but for a larger coupling-strength value, namely β = 10 -3 . It becomes evident that this increase affects the diffusion exponent µ, especially in intervals of K values where AMs are present. In particular, µ gradually decreases to values closer to µ = 1 (normal diffusion) and this decline is more pronounced at larger K values.

Keeping intact the general system's setup (i.e. the kick-strength parameters being equal in all coupled SMs) and using the same ensemble of ICs, we probed the decaying scaling laws of the highest (µ * ) values of the diffusion exponent µ [Eq. ( 14)] (the pronounced peak locations denoted by black filled circles in Fig. 8), as a function of K for increasing β values. These µ * values lie in the first 11 intervals of the kick-strength parameter values K [Eq. ( 4)] for which period p = 1 AMs exist for the 2D SM [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]. The obtained results are shown in Fig. 9. The solid curves in Fig. 9(a) correspond to fittings of the data, obtained after n = 10 4 iterations, with functions µ * = A K B , for each presented case [the considered K values are given in the legend of Fig. 9(d)]. From these curves it becomes evident that strong coupling between neighboring 2D maps tends to evoke a global normal diffusive transport, even when initially the 2D maps may include ensembles of ICs leading to superdiffusive rates due to the presence of AMs. We can also observe that for coupling-strength values β 0.005 the system acquires rather rapidly global normal diffusion rates even for small K values [Fig. 9(a)]. The obtained values of the fitting parameters A and B (and their standard deviations), as a function of the couplingstrength parameter β, are respectively given in Figs. 9(b) and(c). A clear tendency of µ * to become µ * = 1 and independent of K is seen, as A and B respectively become A = 1 and B = 0 [these values are denoted by dashed lines in Figs. 9(b) and (c)] for large enough β values.

Furthermore, and always for the same system setup, we explored the effect of the coupling-strength β, on the highest values µ * of the diffusion exponent µ for these 11 particular K values associated with strong superdiffusion. These results are presented in Fig. 9(d), where it becomes clear that as the couplingstrength β gets stronger the system tends to suppress the superdiffusion and gradually acquires a transport rate µ * ≈ 1. This behavior is typically observed in purely chaotic regions in the case of a single SM and for ICs on grids not containing AMs, for which the diffusion is normal and is characterized by µ ≈ 1

From the analysis of Figs. 8 and9 we observe that the coupling-strength β operates as a control mechanism, being able to suppress superdiffusion rate (µ > 1) into normal diffusion (µ ≈ 1). In Sect. 4.1, we performed an investigation of the relation between the delay to reach ballistic transport (µ = 2) and the kick-strength parameter K (which affects the fraction of chaotic ICs) in single SMs [see Fig. 6(a)]. Here we attempt a similar investigation for the coupled 2ND system (6) trying to associate the number of iterations required for the diffusion exponent µ [Eq. ( 14)] to reach the value µ = 2 with the value of β.

We first plot in Fig. 10(a) the variance N j=1 (∆y j ) 2 [Eq. ( 14)] for the same case of coupled SMs we considered so far in Figs. 8 and9 (i.e. N = 5, fixed K j = K = 6.5) for increasing values of the coupling-strength parameter β. Fig. 10(b) shows the corresponding numerically computed diffusion exponent µ [Eq. ( 14)] as a function of the map's iterations n. One can notice that for all β values the coupled system gradually shifts from superdiffusion (µ > 1) to normal diffusion (µ = 1). This transition happens faster as β grows.

In order to investigate the system's chaoticity we present in

β = 0.00001 β = 0.00005 β = 0.00010 β = 0.00050 β = 0.00100 β = 0.00500 β = 0.01000 β = 0.05000 (c) -16 -12 -8 -4 0 1 2 3 4 log 10 <GALI 2 > log 10 (n)
Figure 10: (a) The variance N j=1 (∆y j ) 2 [Eq. ( 7)] of the coupled SMs system (6) with N = 5 and K j = K = 6.5, j = 1, 2, . . . , 5, for various values of the coupling-strength parameter β (given in the legend) for the same setup and ensemble of ICs described in Figs. 8 and9, as a function of the system's iterations n. The black dotted and dashed lines respectively correspond to µ = 1 and µ = 2. (b) The corresponding numerically computed diffusion exponent µ [Eq. ( 14)] as a function of the map's iterations n. For each β value the coupled system gradually shifts from superdiffusion (µ > 1) to normal diffusion (µ = 1). (c) Average (over ≈ 100, 000 ICs) value Λ of the ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] as a function of β after n = 10 5 iterations. The error bars denote one standard deviation in the computation of the average value. Inset: the evolution of Λ with respect to the number of iterations n. Higher Λ values correspond to larger β, for which normal diffusion rates are also observed. (d) Average GALI 2 values (12) for the same set of ICs as in (c). All GALI 2 decay exponentially fast to zero as the collective dynamics of the particular ensembles of ICs is chaotic. The respective standard deviations (not plotted here) are very small and hardly visible. Fig. 10(c) the average (over all considered ICs) value Λ of the ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] for increasing values of β, after n = 10 5 iterations, while the index's evolution is depicted in the figure's inset. Note that larger Λ values correspond to stronger coupling-strength cases, which exhibit normal diffusion rates. In other words, for this particular setup of the coupled SMs and ensembles of ICs, stronger coupling leads to stronger chaos as measured by the average ftMLE. This is also evident by calculating the average GALI 2 values [Fig. 10(d)], as all GALI 2 exponentially decay to zero, indicating collective chaotic behavior, which is also associated with the positive average ftMLE. Moreover, we observe that this GALI 2 rate of decay is related to the Λ values, i.e. the smaller Λ is the slower GALI 2 tends to zero [compare for example the blue (β = 0.00001) and the brown (β = 0.05) curves in Figs. 10(c) and 10(d)]. We note that the relation between the rate of exponential decay of the GALI 2 and the ftMLE was also mentioned in Sect. 3 (see also [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method[END_REF] for more details).

Trying to further investigate the effect of the couplingstrength parameter β on the chaotic behavior of the coupled SMs we present in Fig. 11 the system's Λ for the K values considered in Fig. 9 and for different β values covering a spectrum of 4 orders of magnitude (from β = 10 -4 to β = 1). We note that Λ is computed as the average ftMLE over 10,000 ICs (i.e. points on a 100 × 100 grid covering the system's entire phase space) after n = 10 5 iterations. From the results of Fig. 11 we see that Λ increases as K grows for fixed β values, following a trend similar to the one seen in Fig. 6(b) for the single 2D SM (denoted by the dashed curve in Fig. 11), having values (in log 10 scale) slightly larger than the ones obtained in that case. Furthermore, a steady increase of log 10 Λ for increasing β values is observed for each considered K value. Thus, for relative large K values in the range 6.5 K 70, increasing the coupling between identical 2D SMs results to chaotic motion described by slightly growing ftMLEs, with respect to the ones observed for the single 2D SM.

Table 2: Nomenclature of different coupled SMs arrangements. In Fig. 13 we present results for the coupled system (6) with N = 5 for which each individual 2D map has a single (S) kick-strength value K j = K = 6.5, j = 1, 2, . . . , 5. Note that K = 6.5 corresponds to the appearance of a p = 1 AM. In all arrangements (named Arrangements A) the considered ICs on the j = 3, 2D map cover the whole phase space (left panels in the two quartets of panels in Fig. 12), while in the other four maps ( j = 1, 2, 4, 5) ICs lie on the phase space regions named 1B, 1C, 1D (see Fig. 12, upper row of Fig. 7 and Table 1) exhibiting different fractions of chaotic motion: respectively P C ≈ 75%, P C ≈ 50% and P C ≈ 25%. In Fig. 14 we show results for a different set of arrangements (named Arrangements B) for which the central j = 3 map corresponds to cases 1B, 1C or 1D and the remaining j = 1, 2, 4, 5 maps contain ICs on the whole phase space. For all the 6 presented cases the coupling-strength parameter is set to β = 0.001. Figure 11: Average [over 10,000 ICs on a 100 × 100 grid covering the system's entire phase space (see text for more details)] value Λ of the ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] of the coupled SMs system [START_REF] Karney | Long-time correlations in the stochastic regime[END_REF] with N = 5 and K j = K, j = 1, 2, . . . , 5, after n = 10 5 iterations, for the K values used in Fig. 9 and for different coupling-strengths β = 10 -4 , 10 -3 , 10 -2 , 10 -1 , 1 (respectively blue, green, orange, purple and red points/curves). The dashed curve corresponds to the law Λ = ln(K/2) seen in the case of a single 2D SM (see Fig. 6(b)].

Case name

Coupled 2D SMs with equal kick-strengths and different fractions of chaos around their AMs

In what follows, we examine the role of the total amount (percentage) of chaotic ICs in the evolved ensembles on the longterm diffusion rate properties for systems of coupled SMs with homogeneous kick-strengths K j = K, j = 1, 2, . . . , 5. The chaoticity and diffusion rate of these particular ensembles of ICs are shown in Fig. 12 for the case of K = 6.5, for which an AM of period p = 1 exists (see Table 1). The first four left panels show phase space color plots based on GALI 2 calculations, while the remaining four panels are color plots based on the value of the diffusion exponent µ [Eq. ( 2)], in a similar fashion to Fig. 7. The first, left, unlabeled panel in each group depicts an ensemble of ICs covering the entire phase space [0, 2π) × [0, 2π), for which extensive chaos is present (P C ≈ 99%). Panels 1B, 1C and 1D refer to ensembles of ICs (the explicit ranges of these regions are given in Table 1) around an AM of period p = 1, containing respectively P C ≈ 75%, P C ≈ 50% and P C ≈ 25% chaotic orbits.

We use these ensembles of ICs to form different coupled SMs (6) Arrangements and we measure again the resulting diffusion exponent µ [Eq. ( 14)] and average chaos indicators Λ and GALI 2 , by setting the coupling-strength parameter at a moderate fixed value, β = 0.001. These choices result to coupled SMs setups which allow us to investigate, in feasible CPU times, the effect of ICs with different chaos percentages on the diffusion properties in the presence of p = 1 AMs in the 2D map components of the coupled systems.

The first setup we considered (named Arrangement A) comprises of coupled 2D SMs with single (S) (i.e. the same) K values, K j = K = 6.5, j = 1, 2, . . . , 5, for which an AM of period p = 1 exists, and the same central ( j = 3) map arrangement for all cases, having ICs covering the entire phase space (see left panels in each quartet of figures in Fig. 12). In the other four maps the ICs lie on regions containing different fractions of chaotic motion, namely P C ≈ 75%, P C ≈ 50% and P C ≈ 25% (these numerical values are used in defining the name of each set up in Table 2). The second setup type (Arrangement B) consists of coupled 2D SMs with again single (S) K values (K j = K = 6.5) for which the j = 3 map corresponds to regions with various fractions of chaotic motion (P C ≈ 75%, P C ≈ 50% and P C ≈ 25%), while the remaining maps are identical having ICs chosen on a grid covering the entire phase space of the 2D map.

Our findings regarding Arrangement A are shown in Fig. 13 where we consider the three cases (referring to different fractions of chaotic motion) described in Table 2. For example, the case with name S75A refers to a system of 5 coupled 2D SMs with single kick-strength values, namely K j = K = 6.5, for which an AM of period p = 1 exists. In this case the ICs in the central ( j = 3) map are on a grid covering the entire phase space, while the other four maps are identical having ICs corresponding to a P C ≈ 75% level of chaotic orbits. Fig. 13(a) shows the evolution of 5 j=1 (∆y j ) 2 [Eq. [START_REF] Mackay | Stochasticity and Transport in Hamiltonian Systems[END_REF]] for ≈ 100, 000 ICs of the coupled system. We can see that only case S25A tends to eventually recover a ballistic diffusion rate, while cases S75A and S50A seem to converge to normal diffusion. In the inset we show the evolution of the diffusion exponent µ [Eq. 14] as a function of the system's iterations n for the same 3 cases. Figure 12: Phase space portraits of the 2D SM ( 1) with K = 6.5 (for which a p = 1 AM exists; see Table 1) colored according to the GALI 2 [START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF] values (left four panels), and the diffusion exponent µ [Eq. ( 2)] (right four panels). The left, unlabeled panel in each quartet of figures depicts an ensemble of ICs on the entire phase space, [0, 2π) × [0, 2π), exhibiting extensive chaos (P C ≈ 99%). Panels labelled as 1B, 1C and 1D (as in Table 1) refer to ensembles of ICs around the AM with respectively P C ≈ 75%, P C ≈ 50% and P C ≈ 25%.

(b)

Figure 13: Cases S75A (blue curves), S50A (green curves) and S25A (red curves) of the coupled SMs system (see Table 2), (a) The variance 5 j=1 (∆y j ) 2 of an ensemble of ≈ 100, 000 ICs, as a function of the system's iterations n. Only case S25A shows a tendency to recover ballistic diffusion rate (as a slope ≈ 2 at the later stages of the evolution indicates), while cases S75A and S50A seem to converge to a normal diffusion rate (slope ≈ 1). Inset: the evolution of the diffusion exponent µ [Eq. ( 14)] as a function of the iterations n for the same cases. (b) The evolution of the average ftMLE Λ as a function of the iterations n for the cases of (a). In all cases Λ converges to approximately the same non-zero value. Inset: the time evolution of the corresponding GALI 2 which decays exponentially to zero as the collective dynamics of these ensembles is chaotic. 2.

The evolution of the average ftMLE Λ as a function of the iterations n is given in Fig. 13(b). We see that all Λ values show a tendency to converge towards a similar, positive value indicating chaotic behavior. An important observation here is that the larger the fraction of chaotic ICs is in the non-central maps the faster the convergence to this asymptotic Λ value takes place. The inset of Fig. 13(b) shows the evolution of the corresponding GALI 2 values which decay exponentially fast to zero as the collective dynamics of this particular ensembles of ICs is chaotic. Figure 15: Phase space portraits of the 2D SM (1) with K = 3.1 (for which a p = 4 AM exists; see Table 1) colored according to the GALI 2 [START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF] values (left four panels), and the diffusion exponent µ [Eq. ( 2)] (right four panels). Panels labelled as 4A, 4B, 4C and 4D (as in Table 1) refer to ensembles of ICs around the AM with respectively P C ≈ 99%, P C ≈ 75%, P C ≈ 50% and P C ≈ 25%.

Table 3: Nomenclature of different coupled SMs arrangements. In Fig. 16 we present results for the coupled system (6) with N = 5, having mixed (M) K values: always K 3 = 3.1 (corresponding to the existence of an AM of period p = 4) for the central map ( j = 3), for which ICs are chosen on a grid covering the 4A region of Fig. 15 (see also Table 1), while the other four maps have K 1 = K 2 = K 4 = K 5 = 6.5 (corresponding to the appearance of a p = 1 AM) with ICs lying on phase space regions 1B, 1C and 1D (see Fig. 12, upper row of Fig. 7 and Table 1) with P C ≈ 75%, P C ≈ 50% and P C ≈ 25% (Arrangements A). In Fig. 17 we show results for another set of arrangements (named Arrangements B) for which the j = 3 map corresponds to a 2D map with K 3 = 6.5 and ICs on its entire phase space (left panels in the two quartets of panels in Fig. 12), and the remaining maps have K 1 = K 2 = K 4 = K 5 = 3.1, belonging to the cases 4B, 4C and 4D (see Fig. 15, lower row of Fig. 7 and Table 1) with P C ≈ 75%, P C ≈ 50% and P C ≈ 25%. For all the 6 presented cases the coupling-strength parameter is set to β = 0.001.

4C 4C [0, 2π) × [0, 2π) 4C 4C M25B 4D 4D (K = 6.5 -AM period: 1) 4D 4D
The investigation of cases belonging to what we call Arrangement B reveals different trends and the corresponding results are presented in Fig. 14. We here consider cases S75B, S50B and S25B (see Table 2 for more information). For instance, case S50B refers to a setup for which again single (S) K values (K = 6.5) are considered with the central ( j = 3) map corresponding to a P C ≈ 50% level of chaos, while the four other maps have ICs on the entire phase space (see also Fig. 12). Fig. 14(a) depicts results similar to the ones of Fig. 13(a). All cases show a common tendency to converge to normal diffusion rate (µ ≈ 1), something which is also seen in the inset of Fig. 14(a), as the diffusion exponent µ of these particular cases does not show signs of approaching an asymptotic extreme rate (ballistic transport), at least within the considered time scales. Regarding the average global chaoticity of these ensembles of ICs, we observe a more uniform evolution both for the ftMLEs Λ [Fig. 14(b)], which converge to approximately the same non-zero value, and the GALI 2 values [inset of Fig. 14(b)], which follow very similar exponential decay rates.

Coupled 2D SMs with different kick-strengths and AMs' periods

Let us now extend our analysis by considering alternative coupled SMs configurations, namely setups with different kickstrength values for individual 2D maps containing AMs of different periods. Similarly to the study performed in Sect. 4.2.2, we explore diffusion rate properties and the system's global chaoticity for some cases belonging to what we name Arrangements A and B, including SMs with different fractions of chaotic ICs and various [mixed (M)] values of kick-strengths. The de-tails of these configurations are given in Table 3. Arrangements A contain in the off-center maps ( j = 1, 2, 4, 5) ensembles on regions 1B, 1C and 1D (see Fig. 12, upper row of Fig. 7 and Table 1), while arrangements B include the cases 4B, 4C and 4D presented in the upper row of Fig. 7 and in Table 1. In addition, color plots of µ and GALI 2 values for the latter 2D map regions are seen in Fig. 15. More specifically, region 4A corresponds to an ensemble of ICs within the intervals [0.0, 1.047] × [1.964, 3.142], where an AM of period p = 4 exists, and exhibits extensive chaos (P C ≈ 99%). Regions 4B, 4C and 4D refer to ensembles of ICs around the same AM of period p = 4, respectively containing P c ≈ 75%, P c ≈ 50% and P c ≈ 25% of chaotic ICs.

The results for Arrangement A (where the central map contains an AM of period p = 4, while the other maps contain AMs of period p = 1) are presented in Fig. 16. From the results of Fig. 16(a) we see that all cases initially exhibit superdiffusive spreading (µ ≈ 2), while after some iterations (n ≈ 10 5 ) they show a tendency to halt this process, exhibiting very low subdiffusion rates (µ → 0). However, towards the end of the simulation the respective rates show a tendency to increase again. These changes are better depicted in the inset of Fig. 16(a) where the evolution of the diffusion exponent µ [Eq. [START_REF] Kroetz | Hidden high period accelerator modes in a bouncer model[END_REF]] is given as a function of the number n of iterations. The respective ftMLE Λ converges to approximately the same positive value for all cases [Fig. 16(b)]. Analogous trends for Λ to the ones presented in Fig. 13(b) are also observed here: the larger the fraction of chaotic ICs is in the non-central maps, the faster the convergence to the index's asymptotic value happens. In

(b)

Figure 16: Similar to Fig. 13 but for cases M75A (blue curves), M50A (green curves) and M25A (red curves) of Table 3. 15, lower row of Fig. 7 and Table 1), while the ICs in the central map with K 3 = 6.5 cover the whole phase space. For these configurations, we find that all cases initially exhibit superdiffusion rates (µ ≈ 1.5) [Fig. 17(a) and its inset], which are lower than the ones seen in Arrangement A cases [Fig. 16(a)], while after some iterations (n ≈ 10 5 ) they show a tendency to gradually decrease these rates towards normal diffusion behaviors (µ → 1). Fig. 17(b) renders the respective time evolution of the average ftMLE Λ for these cases. Again all Λ values seem to converge to the same asymptotic limit with the fraction P C of chaotic ICs of the central map practically not affecting the route to this limiting value. The dynamical proximity of all these cases is also seen in the evolution of GALI 2 in the inset of Fig. 17

Summary and discussion

In this work, we examined the long-term diffusion and chaoticity properties of single 2D SMs, as well as systems of 5 coupled 2D maps, focusing our attention on parameter values for which the phase space of the respective systems exhibit anomalous diffusion rates due to the presence of AMs of different periods.

For single 2D SMs we reviewed the most typical diffusion properties for regular (subdiffusion) and chaotic (normal diffusion) orbits, as well as stable/unstable AMs of different periods (with the former exhibiting ballistic transport). For the considered sets of ICs, we have also quantified their chaoticity, using the GALI 2 [START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF] and ftMLE [START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] indices. Furthermore, we investigated the map's diffusion properties, by measuring the diffusion exponent µ [Eq. [START_REF] Rechester | Calculation of turbulent diffusion for the Chirikov-Taylor model[END_REF]] and the effective diffusion coefficient D eff [START_REF] Meiss | Correlations of periodic, area-preserving maps[END_REF], for a range of kick-strength values K and for different final iteration numbers (Fig. 3). A systematic study of the impact of the number of iterations, in the presence of AMs of period p = 1, lead to the conclusion that all AMs of period p = 1 found in the interval K ∈ [0, 70] asymptotically acquire an extreme diffusion rate (ballistic transport with µ ≈ 2), even when for relatively short times they behave closely to normal diffusion (µ ≈ 1, see Fig. 5). We have found that the larger the K value of the SM with a p = 1 AM is, the more chaotic is its respective phase space as measured by the average ftMLE Λ for large ensembles of ICs [Fig. 6(b)], with Λ being accurately described by the law Λ = ln(K/2) [START_REF] Shevchenko | On the maximum Lyapunov exponent of the motion in a chaotic layer[END_REF][START_REF] Harsoula | Characteristic times in the standard map[END_REF]. Moreover, we considered SMs with kick-strength values where AM of different periods occur. For these cases, we investigated the impact of the choice of the ensemble of ICs (i.e. practically changing the fraction P C of chaotic ICs around the AM) on the diffusion exponent µ. We found that the more chaotic the ensemble is, the longer it takes for ballistic transport to occur, i.e. for the diffusion exponent µ to converge to µ = 2 (Fig. 7).

In the second part of our analysis, we extended our investigation by exploring diffusion and chaos properties for coupled SMs containing AMs of different periods. We began our analysis by choosing identical ensembles of ICs on a grid covering the entire phase space of each coupled 2D map, with equal kick-strength values K and moderate couplings (β = 10 -4 and 10 -3 ) for all maps, and identified the main K intervals where AMs of different periods appear (Fig. 8). We derived scaling laws for the regions where AMs of period p = 1 occur (pronounced peaks in Fig. 8), which describe the dependence of those peaks on K [Figs. 9(a)-(c)], and estimated the respective global diffusion exponent µ * for a range of coupling parameter values β [Fig. 9(d)]. For AMs of period p = 1, the main finding is that, as β increases the global diffusion exponent decreases, tending towards µ * = 1 with a faster rate as K reaches larger values. By setting the kick-strength values at K j = K = 6.5, j = 1, 2, . . . , 5 for all SMs (corresponding to the occurrence of the strongest superdiffusion in each 2D map, with µ ≈ 2), we obtained a global asymptotic diffusion for the coupled system characterized by normal diffusion rates (µ ≈ 1) [Fig. 10(b)]. Furthermore, this convergence towards normal diffusion rates takes place faster as the coupling parameter β increases. The trend of the relatively faster convergence to normal diffusion is also found to be related to the global chaoticity (as measured by the ftMLE Λ and the GALI 2 ) of the ensembles of ICs. Namely, the average ftMLE value increases (indicating stronger chaos) monotonically with the increase of the coupling parameter β [Figs. 10(c) and(d)].

We have also performed a similar analysis for different arrangements of the 5 coupled SMs system, in order to investigate how different configurations of individual 2D maps affect the respective global diffusion rates (see Tables 2, and3). Namely, we sought to find conditions under which the initially ballistic transport (due to the presence of AM of low periods) can be suppressed by the presence of neighboring maps without AMs and ensembles of chaotic ICs exhibiting normal diffusion rates. In order to achieve this goal we performed extensive numerical simulations for two SMs' arrangements.

In the Arrangement A type of coupled maps we considered single (presence of AMs of period p = 1 in all SMs, see Fig. 13) and mixed (presence of AMs of period p = 4 for the central, j = 3, SM and of period p = 1 for the others, see Fig. 16) K values for fixed coupling-strength β = 0.001. In both cases the diffusion exponent starts at µ ≈ 2 (due to the presence of the AMs), then reaches a plateau characterized by µ ≈ 0, and finally recovers the initial rate in the single case for ensembles with relatively low fraction P C of chaotic ICs (case S25A), while it tends to normal diffusion values for relatively high P C fractions (cases S75A and S50A). For the mixed cases (M75A, M50A, M25A) this trend is not clearly evident [Fig. 16(a)] for the maximum number of iterations we managed to reach. In terms of global chaoticity (quantified through ftMLE Λ computations), we found that all Λ values show a tendency to converge towards a similar positive value. The larger the fraction P C of chaotic ICs is in the off-center maps, i.e. respectively cases S75A, S50A, S25A and M75A, M50A, M25A the faster the convergence to this value takes place.

In the Arrangement B type of coupled maps we considered single (central SM with an AM of period p = 1 and varying fractions of chaotic ICs, while the remaining SMs have ICs in the whole phase space, see Fig. 14) and mixed (central SMs with a p = 1 AM and ICs in the whole phase space, along with varying fractions P C for the other SMs, which have AMs of period p = 4, see Fig. 17) K values with β = 0.001. For both cases we found that the diffusion rate starts again at being strongly anomalous with µ ≈ 2 (due to the presence of AMs of period p = 1), while it gradually converges to a normal diffusion rate with µ ≈ 1 (cases S75B, S50B, S25B and M75B, M50B, M25B). Furthermore, in all cases, the average ftMLE Λ asymptotically tend to similar values [see Fig. In summary, the only configuration we found (up to the maximum number of iterations we managed to reach) in which the system globally recovers extreme diffusion rates is case S25A in Fig. 13(a). In all other explored cases, we observed a suppression of the diffusion rate. Thus, this case deserves further investigation for larger number of iterations, in order to determine possible changes in the diffusion process. This is a task we plan to tackle in the future.

Let us also mention here that in our analysis we did not investigate coupling combinations of SMs that may include weakly chaotic or sticky ICs, which may be confined or bounded by invariant tori, and whose diffusion can be more complex (potentially exhibiting subdiffusion rates). Instead, we focused our attention only on cases where chaotic motion dominates the phase space, i.e. the chosen 2D SMs have relatively large kickstrength values. Taking into account SMs and ensembles of ICs leading to subdiffusion and coupling them together will most likely result in a much richer global diffusion behavior, which deserves separate investigations. Finally, we note that we limited our study to coupled SMs with a relatively small number, N = 5, of 2D maps. This choice allowed us to analyze rather large ensembles of ICs, as well as to perform extensive numerical simulations and calculate diffusion exponents and chaos indicators for large numbers of map iterations. facility of the University of Cape Town, as well as the Centre for High Performance Computing (CHPC) of South Africa for providing the computational resources needed for obtaining the numerical results of this work.

Figure 1 :

 1 Figure 1: Results for the 2D SM (1) with [(a), (b) and (c)] K = 3.1 and [(d), (e) and (f)] K = 6.5. In (a) and (d) phase space portraits of the SM are shown, which are produced by 10 6 iterations of several ICs (black points), while colored points are used to identify particular orbits. In (a) the consequents of the chaotic orbit with IC (x 0 , y 0 ) = (2, 1.5) (red pentagon point) are plotted in red, while blue points correspond to a regular orbit with IC (x 0 , y 0 ) = (3, 0.5) (blue diamond point). In (d) the location of two AMs of period p = 1 are identified, a stable AM (green circle point) at (x 0 , y 0 ) = (1.8298, 0) and an unstable AM (purple square point) at (x 0 , y 0 ) = (1.3118, 0). The evolution of the GALI 2 of all these orbits, with respect to the number of iterations n, is illustrated in (b) for the chaotic (red curve) and the regular (blue curve) orbit of (a), and in (e) for the stable (green curve) and the unstable AM (purple curve) of (d). In (b) GALI 2 (n) tends to zero exponentially fast for the chaotic orbit, while it decreases to zero following a power law GALI 2 ∝ n -2 (indicated by a dashed line) for the regular orbit. In (e) GALI 2 (n) remains practically constant for the stable AM, while it tends to zero exponentially fast for the unstable AM. In (c) and (f) we see the evolution of the ftMLE Λ(n)[START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF], respectively for the chaotic (red curve) and the regular (blue curve) orbit of (a), and for the stable (green curve) and the unstable (purple curve) AM of (d). Λ(n) saturates to a positive value in (c) and (f), respectively for the chaotic orbit and the unstable AM, while it tends to zero following a power law Λ(n) ∝ n -1 (indicated by a dotted line in both panels) for the regular orbit and the stable AM.

Figure 2 :

 2 Figure2: (a) The variance (∆y)2 (2) of the angular momentum y n of the 2D SM (1) as a function of the map's iterations n, for orbits in four different small regions of the map's phase space. For each region a set of ICs on a 315 × 315 grid of equally spaced mesh points (i.e. in total ≈ 100, 000 ICs) are taken and the variance is obtained as average over all these orbits. For K = 3.1 ICs in the region [2.995, 3.005] × [0.495, 0.505] (i.e. 2.995 ≤ x 0 ≤ 3.005 and 0.495 ≤ y 0 ≤ 0.505) and [1.995, 2.005] × [1.495, 1.505] containing respectively the IC of the regular (blue curve) and the chaotic (red curve) orbit of Fig.1(a)are considered, while the green and purple curves respectively correspond to regions[1.825, 1.835] × [0.0, 0.01] [containing the stable AM of Fig.1(d)] and[1.295, 1.315] × [0.0, 0.01] [surrounding the unstable AM of Fig.1(d)] for the K = 6.5 SM. The dashed and dotted lines respectively indicate functions proportional to n 2 and n 1 . Note that the red and purple curves overlap. (b) The numerically computed diffusion exponent µ of (2) for the curves of (a), where the same color for each case is used.

  Figure2: (a) The variance (∆y)2 (2) of the angular momentum y n of the 2D SM (1) as a function of the map's iterations n, for orbits in four different small regions of the map's phase space. For each region a set of ICs on a 315 × 315 grid of equally spaced mesh points (i.e. in total ≈ 100, 000 ICs) are taken and the variance is obtained as average over all these orbits. For K = 3.1 ICs in the region [2.995, 3.005] × [0.495, 0.505] (i.e. 2.995 ≤ x 0 ≤ 3.005 and 0.495 ≤ y 0 ≤ 0.505) and [1.995, 2.005] × [1.495, 1.505] containing respectively the IC of the regular (blue curve) and the chaotic (red curve) orbit of Fig.1(a)are considered, while the green and purple curves respectively correspond to regions[1.825, 1.835] × [0.0, 0.01] [containing the stable AM of Fig.1(d)] and[1.295, 1.315] × [0.0, 0.01] [surrounding the unstable AM of Fig.1(d)] for the K = 6.5 SM. The dashed and dotted lines respectively indicate functions proportional to n 2 and n 1 . Note that the red and purple curves overlap. (b) The numerically computed diffusion exponent µ of (2) for the curves of (a), where the same color for each case is used.

Figure 3 :

 3 Figure 3: Diffusion exponent µ of (2) [(a) and (c)] and effective diffusion coefficient D eff (5) [(b) and (d)] as a function of K after n = 10 4 [(a) and (b)] and n = 10 5 [(c) and (d)] iterations of the 2D SM (1) for a set of 315 × 315 ICs on a grid of equally spaced mesh points (≈ 100, 000 ICs) on the map's entire phase space [0, 2π) × [0, 2π). The intervals on the black horizontal line at the bottom of (a) and (c) denote the existence intervals [Eq. (4)] of period p = 1 AMs. The high peaks in (a) and (c) for K > 2π (their highest values µ * are marked by black circles) correspond to cases of anomalous diffusion due to the presence of period p = 1 AMs, while regions with µ ≈ 1 tally to normal diffusion processes. Intervals of K values related to anomalous diffusion also generate the peaks in (b) and (d).

Figure 4 :

 4 Figure4: The variance (∆y)2 (a) and the corresponding numerically computed diffusion exponent µ [Eq. (2)] (b) of the 2D SM (1) as a function of the map's iterations n, for an ensemble of orbits with ICs on a 315 × 315 grid covering the map's entire phase space [0, 2π) × [0, 2π), and four K values, namely K = 6.866 (blue curves), K = 25.291 (green curves), K = 44.002 (red curves) and K = 69.173 (black curves), which correspond to some of the black circle denoted peaks in Fig.3(a). The black dotted and dashed lines in (a) respectively correspond to µ = 1 and µ = 2.

Figure 5 :

 5 Figure 5: (a) The highest values µ * of the diffusion exponent µ [Eq. (2)] in the first 11 intervals of K values [Eq. (4)] for which period p = 1 AMs exist for the 2D SM (1), as a function of K for n = 10 3 [blue diamonds], n = 10 4 [green pentagons; points denoted by black circles in Fig. 3(a)], n = 10 5 [red squares; points denoted by black circles in Fig. 3(c)], n = 10 6 [orange triangles] and n = 10 7 [purple circles] iterations of the map. In all cases ensembles of orbits with ICs on a 315 × 315 grid covering the map's entire phase space [0, 2π) × [0, 2π) were considered. The solid curves correspond to fittings of the data with functions µ * = AK B for each presented case. The obtained values of the fitting parameters A and B (along with their uncertainties), as a function of the number of iterations n are respectively given in (b) and (c). The values A = 2 in (b) and B = 0 in (c) are indicated by horizontal dashed lines.

Figure 6 :

 6 Figure 6: (a) The number n of iterations of the 2D SM (1) needed for the diffusion exponent µ [Eq. (2)] to reach the µ = 2 value (ballistic transport) for K values at which µ becomes maximum [see Figs. 3(a), (c) and Fig. 5(a)] in the first 11 intervals of existence [Eq. (4)] of period p = 1 AMs, as a function of K (and in the inset as a function of the width, ∆K l = (2πl) 2 + 16 -2πl, of these intervals). (b) Average value Λ [over 10,000 orbits with uniformly distributed ICs on a 100 × 100 grid covering the map's entire phase space [0, 2π) × [0, 2π)] of the orbits' ftMLE[START_REF] Ouchi | Advective Diffusion of Particles in Rayleigh-Bénard Convection[END_REF] after n = 10 5 iterations, for the same values of K (explicitly given in the inset) as in (a) The error bars denote one standard deviation in the computation of the average value, while the dashed red curve corresponds to the law Λ = ln(K/2). Inset: the evolution of Λ with respect to the map's iteration n. Higher Λ values correspond to larger K values.

KKKFigure 7 :

 7 Figure 7: Phase space portraits of the 2D SM (1) colored according to the GALI 2[START_REF] Mori | Energy Dissipation and Its Fluctuations in Chaotic Dynamical Systems[END_REF] value (panels in the two left columns) after n = 50 iterations of orbits with ICs on a dense grid covering the considered phase space region, and the diffusion exponent µ [Eq. (2)] (panels in the two middle columns) computed by appropriately fitting the variance (∆y)2 (2) values for iterations n in the interval 10 3 ≤ n ≤ 10 4 (see text for more details). Each row of panels corresponds to a different K value, namely K = 6.5, K = 4.0844, K = 6.9115 and K = 3.1 (from top to bottom), for which stable AMs of period p = 1, 2, 3 and 4 respectively exist. Points are colored according to the color scales on top of the panels. In the GALI 2 plots yellow regions refer to regular motion, dark blue and black areas correspond to chaotic motion, while in-between colors to weakly chaotic ICs requiring more iteration to reveal their chaotic nature. In the µ color plots ICs leading to subdiffusion rates are colored in dark blue (0 ≤ µ < 1), normal diffusion cases are depicted in dark red (µ ≈ 1), superdiffusion corresponds to lighter shades of red (1 ≤ µ < 2) and ballistic transport is denoted by yellow (µ ≈ 2). The positions of the stable p = 1 [(1.8298, 0)], p = 2 [(4.211, 4.9324)], p = 3 [(1.80, 0.0)] and p = 4 [(0.3486, 2.6121)] AMs are denoted by small black rectangular points in the right GALI 2 and µ panels. In each case 6 different regions of the phase space are considered. The names of these regions, along with their ranges in the x and y coordinates, are given in Table1, and their positions are indicated by different rectangles in the color plots accompanied by the names of the regions. These regions are ordered in Table1in decreasing content of chaotic orbits (quantified by the percentage P C of ICs leading to chaotic motion) and they were chosen so that in all cases ensembles with similar names have equivalent P C values, i.e. P C ≈ 99% for 'A' regions, P C ≈ 75% for 'B', P C ≈ 50% for 'C', P C ≈ 25% for 'D', P C ≈ 12.5% for 'E' and P C ≈ 6.5% for 'F'. The panels in the last column depict the diffusion exponent µ [Eq. (2)] of the 6 different regions as a function of the map's iterations n.

Figure 8 :

 8 Figure 8: Similar to Fig. 3 but for the N = 5 coupled system (6) with identical kick-strength values K j = K, j = 1, 2, . . . , 5. (a) Diffusion exponent µ [Eq. (14)] as a function of K for fixed coupling-strength β = 10 -4 after n = 10 4 iterations. The considered ensemble of orbits consists of ICs whose projections in each one of the 5 2D SMs are on a 315 × 315 equally spaced grid covering the whole map (see text for more details). In the inset, we plot the effective diffusion coefficient D N eff (7) as a function of K. (b) Same as panel (a) but for a stronger coupling, β = 10 -3 .

Figure 9 :

 9 Figure 9: (a) The highest values µ * of the diffusion exponent µ [Eq.[START_REF] Kroetz | Hidden high period accelerator modes in a bouncer model[END_REF]] in the first 11 intervals (see also Fig.8) of K values [Eq. (4)] for which period p = 1 AMs exist for the 2D SM (1), as a function of K for various coupling-strength β values (see legend), computed for n = 10 4 iterations of the multidimensional map[START_REF] Karney | Long-time correlations in the stochastic regime[END_REF] with N = 5 (see text for more details on the used ensemble of ICs). The solid curves correspond to fittings of the data with functions µ * = A K B for each presented case. The obtained values of the fitting parameters A and B (along with their uncertainties), as a function of β are respectively given in (b) and (c). The values A = 1 in (b) and B = 0 in (c) are indicated by horizontal dashed lines. (d) The highest values µ * of the diffusion exponent µ as a function of the coupling-strength parameter β for different (but fixed in all maps) values of the kick-strength parameter K j = K, j = 1, 2, . . . , 5, for which pronounced superdiffusion transport takes place (corresponding to the K values indicated by black filled circles in Fig.8).

Figure 14 :

 14 Figure 14: Similar to Fig. 13 but for cases S75B (blue curves), S50B (green curves) and S25B (red curves) of Table2.

  M50B

Figure 17 :

 17 Figure 17: Similar to Fig. 13 but for cases M75B (blue curves), M50B (green curves) and M25B (red curves) of Table3

  (b).

  14(b) and Fig.17

  (b)], but in a more uniform manner than in the cases of arrangement A.

  , namely K = 6.866, 12.877, 19.059, 25.291, 31.543, 37.805, 44.002, 50.345, 56.619, 62.895 and 69.173. In particular, we evolve 10, 000 orbits with ICs on a 100 × 100 grid on the entire phase

6. Acknowledgements H. T. M. acknowledges support by a PhD Fellowship from the Science Faculty of the University of Cape Town and partial funding by the UCT Incoming International Student Award. H. T. M. and Ch. S. thank the High Performance Computing 7. Author contributions H. T. M.: Software, Formal analysis, Investigation, Writing -Original draft, Visualization. T. M.: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing -Original draft, Reviewing and Editing. Ch. S.: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing -Reviewing and Editing, Supervision, Project administration, Funding acquisition.