Power-law trapping in the volume-preserving Arnold-Beltrami-Childress map

Understanding stickiness and power-law behavior of Poincaré recurrence statistics is an open problem for higher-dimensional systems, in contrast to the well-understood case of systems with two degrees-of-freedom. We study such intermittent behavior of chaotic orbits in three-dimensional volume-preserving systems using the example of the Arnold-Beltrami-Childress map. The map has a mixed phase space with a cylindrical regular region surrounded by a chaotic sea for the considered parameters. We observe a characteristic overall power-law decay of the cumulative Poincaré recurrence statistics with significant oscillations superimposed. This slow decay is caused by orbits which spend long times close to the surface of the regular region. Representing such long-trapped orbits in frequency space shows clear signatures of partial barriers and reveals that coupled resonances play an essential role. Using a small number of the most relevant resonances allows for classifying long-trapped orbits. From this the Poincaré recurrence statistics can be divided into different exponentially decaying contributions which very accurately explains the overall power-law behavior including the oscillations.

Understanding stickiness and power-law behavior of Poincaré recurrence statistics is an open problem for higher-dimensional systems, in contrast to the well-understood case of systems with two degrees-of-freedom. We study such intermittent behavior of chaotic orbits in three-dimensional volume-preserving systems using the example of the Arnold-Beltrami-Childress map. The map has a mixed phase space with a cylindrical regular region surrounded by a chaotic sea for the considered parameters. We observe a characteristic overall power-law decay of the cumulative Poincaré recurrence statistics with significant oscillations superimposed. This slow decay is caused by orbits which spend long times close to the surface of the regular region. Representing such long-trapped orbits in frequency space shows clear signatures of partial barriers and reveals that coupled resonances play an essential role. Using a small number of the most relevant resonances allows for classifying long-trapped orbits. From this the Poincaré recurrence statistics can be divided into different exponentially decaying contributions which very accurately explains the overall power-law behavior including the oscillations.

I. INTRODUCTION
Hamiltonian systems generically have a mixed phase space in which regular and chaotic motion coexist. A typical chaotic trajectory is often trapped intermittently close to regular regions of phase space for arbitrary long but finite times. The existence of such sticky behavior is known to have important consequences for chaotic transport, e.g. the existence of anomalous kinetics, Lévy processes, and Lévy flights and finds interesting applications in a variety of physical contexts such as chaotic advection in fluids [1], in the three-body problem [2], driven coupled Morse oscillators [3], confinement in particle accelerators [4], and plasma physics [5], to name a few. The behavior of stickiness in a mixed phase-space can be characterized by the cumulative Poincaré recurrence statistics P (t), which is the probability that a chaotic orbit has not returned to an initial region until time t [6]. Fully chaotic systems typically show an exponential decay of P (t) [7][8][9][10]. However, systems with a mixed phase space usually display a much slower decay, usually well-described by a power-law [2,. This so-called power-law trapping has striking implications for transport in many physical systems [2,3,[36][37][38][39].
Power-law trapping is well-understood for Hamiltonian systems and maps with two degrees-of-freedom [15,17,32,[40][41][42][43][44][45][46]. For example in area-preserving maps, regular tori are one-dimensional and therefore act as complete barriers to transport in phase space. If such tori break up they may turn into so-called cantori which provide partial barriers with slow transport across them. Usually there is a whole hierarchy of such partial barriers which is the origin of the power-law trapping. In contrast, in higher-dimensional systems regular tori have an insufficient dimension to provide complete barriers to transport. For example for a 4d symplectic map, the regular tori are two-dimensional and therefore cannot divide the phase space into dynamically disconnected regions.
Thus already purely for topological reasons additional types of transport are possible, leading for example to the famous Arnold diffusion [47][48][49].
Still, power-law trapping is observed in higherdimensional systems [20,28,29,33,35], but not due to a generalized hierarchy [33,35], so that the mechanism of stickiness is different from that in two-dimensional areapreserving maps. While some generalizations have been proposed [29,[50][51][52], a profound understanding of powerlaw trapping in higher-dimensional systems remains a significant open question [46]. This motivates to investigate power-law trapping for 3d volume-preserving maps, whether features of 2d maps still apply, while those present for 4d maps might already become relevant. For example in a 3d volume-preserving map there can be 2d regular tori which therefore are characterized by two frequencies, as in 4d symplectic maps. On the other hand, the 2d tori can provide absolute barriers to the motion in a 3d volume-preserving map as in 2d area-preserving maps. The dynamical properties of 3d volume-preserving maps have been studied in much detail, see e.g. [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69]. However, there are only a few studies of power-law trapping in 3d volume-preserving maps, see Refs. [70][71][72]. In Ref. [70] it is argued that stickiness in their model originates around the hyperbolic invariant sets embedded in the sticky domain. In Ref. [71] recurrence times statistics of an extended standard map have been analyzed. Power-law decay in the presence of accelerator modes is found in Ref. [72].
In this paper, we investigate power-law trapping in three-dimensional volume-preserving maps using the specific example of the Arnold-Beltrami-Childress (ABC) map. The ABC map is a discretized version of the ABC flows with time-periodic forcing [73,74]. The ABC flows are a class of non-turbulent flows which appear as a solution of the three-dimensional Euler equation and is used as a simple model for the study of chaos in threedimensional steady incompressible flows. The 3d phase space of the flow and the corresponding map is known to have invariant surfaces and chaotic streamlines. The flow was first introduced by Arnold [75] and later Childress [76] demonstrated that these are important models for the kinematic dynamo effect in astrophysical plasma, i.e., a magnetic field can be generated and sustained by the motion of an electrically conducting fluid. The flow also serves as a prototypical model of a force-free magnetic field to study spatial diffusion of charged particles [77]. The ABC map exhibits most of the basic features of a typical three-dimensional, time-periodic, volume-preserving ABC flow. The map has a mixed phase space with periodic, quasi-periodic, and chaotic behavior in the 3d phase space. For specific parameters of the map, the phenomenon of resonance-induced diffusion has been observed which leads to global transport in phase space [74]. The map has been employed as a representative of volume-preserving flows to study chaotic advection of inertial particles in bailout embedding models [78][79][80].
The outline of the paper is as follows. Section II describes the ABC map, its basic properties, and illustrates regular orbits in the 3d phase space in Sec. II A. The variation of the regular region is illustrated in Sec. II B using a sequence of 2d plots of the finite-time Lyapunov exponent. The frequency space of the system is introduced in Sec. II C and a boundary cylinder, separating regular and chaotic motion, is identified. In Sec. III, the stickiness near the regular region is studied using the Poincaré recurrence statistics and power-law trapping has been explained in Sec. III A. As illustration, one longtrapped orbit is displayed in the 3d phase space, on a 2d Poincaré section, and in frequency-time representation in Sec. III B. A quantitative analysis of the approach to the boundary cylinder is done in Sec. III C, demonstrating the significant role played by coupled resonances. Furthermore, in Sec. III D we use the frequency analysis to identify a small number of the most relevant resonances. This allows for classifying long-trapped orbits and by this to divide the Poincaré recurrence statistics into different exponentially decaying contributions. This very accurately explains the overall power-law behavior including the oscillations. In Sec. III E possible partial barriers are briefly discussed. Section IV gives a summary and outlook.

II. ARNOLD-BELTRAMI-CHILDRESS MAP
The ABC map [74] is the discrete-time dynamical system (x n , y n , z n ) → (x n+1 , y n+1 , z n+1 ) given by where x n , y n , z n ∈ [0, 2π[ and periodic boundary conditions are imposed. The map is volume preserving for any choice of the real parameters A, B, and C. The dynamics of the system can be very different, showing a mixed phase space or chaotic motion [55,74]. For the rest of the paper, we fix the parameters (A, B) = (2.0, 1.5) and consider different values of C.
For C = 0 the dynamics in (x, z) becomes independent of the motion in the y-coordinate and therefore reduces to a 2d area-preserving subsystem, (2) Figure 1 shows the 2d phase space for this map for (A, B) = (2.0, 1.5), which has the typical structure of a mixed phase space with regular and chaotic motion.
The fixed point at (x, z) = (3π/2, π) is elliptic and therefore surrounded by regular orbits forming invariant tori (full curves). Embedded within this regular region one has a prominent period-4 regular island. In addition to this central regular island one has a symmetry-related island around the elliptic fixed point at (x, z) = (π/2, 0). Outside of these regular islands one has a large area with predominantly chaotic motion, as illustrated by a longer orbit in the figure. Due to the symmetry, dynamical properties, such as stickiness at the islands are identical, so that we can restrict to the central region around the fixed point at (x, z) = (3π/2, π).

A. 3D phase space
When considering the full 3d map (1) for C = 0, one has in addition to the 2d dynamics (2) the motion in y-direction given by y n+1 = y n + B sin(x n+1 ) + A cos(z n ).
As the 2d system (2) decouples from the motion in ydirection when C = 0, the stable fixed points and periodic points of the 2d map turn into elliptic 1d invari-ant lines for the 3d map. The surrounding regular tori turn into 2d regular tori forming concentrically nested cylinders. This is illustrated in Fig. 2(a). Unstable fixed points and unstable periodic points of the 2d map turn into hyperbolic 1d invariant lines. For C > 0 the dynamics in the (x, z) subsystem couples with the motion in y-direction. Starting from C = 0 and increasing C, the lines of the 1d tori begin to distort gradually and the surrounding cylinders of 2d tori also deform. This is illustrated in Fig. 2(b) for the 2d tori and in Fig. 3 where five 1d tori are shown, which for C = 0 form straight lines in the 3d phase space and become distorted with increasing C. We obtain these 1d tori by the "contraction method" starting from 2d tori and iteratively reducing the radius in the (x, z) plane [81].
Such elliptic 1d tori are surrounded by 2d tori. If these form a cylinder extending along the whole y-direction, they act as a full barrier to the motion, i.e. orbits started inside will never be able to move outside, and vice versa. This is therefore similar to the case of 2d area-preserving maps, where 1d tori form absolute barriers to the motion [42,46]. With increasing perturbation, more and more of these 2d tori break up and resonances can lead to gaps in the families of 2d tori, see Sec. II C below.
A representation similar to Fig. 1 can be obtained for the 3d map by using a 2d slice at y = π: only those points (x n , z n ) of an orbit for which |y n − π| < ε = 10 −3 (4) holds, are plotted, see Fig. 4 for C = 0.06. The parameter ε determines the resolution of the resulting plot and smaller values require more iterations to obtain the same number of points of an orbit in the slice. While there are orbits which never fulfill the slice condition (4), a lot of the geometry of the regular orbits can be seen in this plot. Usually 1d tori correspond to single points in the 2d slice and 2d tori lead to one or several loops. Thus the visual appearance is similar to that of the 2d map shown in Fig. 1. In particular there still seems to be some outer regular curve which corresponds to a 2d cylinder in the 3d phase space. Numerically no apparent holes in this surface have been found. Thus such an outer 2d torus forms an absolute barrier separating the regular region from the chaotic sea. Therefore it can be considered as an analogue of the boundary circle in 2d maps [82].
We will refer to this as boundary cylinder in the following, see Sec. III for further discussion on its numerical determination. Additionally, inside the region enclosed by the boundary cylinder there are some chaotic regions, of which the most prominent one is shown as green dots. It seems that none of these inner chaotic regions is dynamically connected to the outer chaotic sea. Several smaller regular islands are embedded within the chaotic region, some of them are shown in Fig. 4. The orbit in blue is a chaotic orbit which stays outside of the region enclosed by the boundary cylinder.

B. Finite-time Lyapunov exponent
To get an overview of how the regular structures are affected when increasing the coupling C we use the finitetime-Lyapunov exponent (FTLE) as chaos indicator [83][84][85]. The FTLE estimates the rate of separation of orbits started at nearby initial conditions. Large values of the FTLE indicate chaotic motion. Figure 5 shows the FTLE for a sequence of values for C, computed in the plane y = π and a grid of 1000 × 1000 points for x ∈ [π, 2π] and z ∈ [π/2, 3π/2], corresponding to the surrounding of the central regular region. For C = 0, regular motion with small FTLE appears as a large island. Already for C = 0.02 the region with small FTLE has substantially reduced and a significant outer part of the former regular region has turned into a region with intermediate values of the FTLE. Typically, initial conditions in this region will stay there for long time before they eventually leave towards the chaotic sea. This sticky region is of particular This gradual transition is also reflected in the sequence of histograms of the FTLE shown in Fig. 6. For C = 0.02 one has a large peak at small values of the FTLE indicating the large regular region. Chaotic motion corresponds to the peak at larger values and the peak in between corresponds to the intermediate type of motion surrounding the central regular part seen in Fig. 5(b). With increasing C this peak moves towards larger values of the FTLE, indicating more chaotic motion, and eventually merges at C = 0.1. with the peak associated with the chaotic sea. In addition, the fraction of regular orbits substantially decreases with increasing C.
Interestingly, the way how the 2d tori break up for the considered ABC map leads to a quite large chaotic layer inside the regular island, see Fig. 4, which is responsible for the additional peak in the FTLE histograms. This seems to be more pronounced than that for 1d tori in a typical area-preserving map. For example the FTLE distribution for the standard map makes a transition from bimodal to unimodal (Gaussian) when increasing the chaoticity in the system [84]. Also for other parameters the histograms appear essentially bimodal.

C. Frequency space
The frequency-space representation is an important complementary approach to investigate the organization of regular motion in higher-dimensional dynamical systems. With every regular 2d torus in the 3d phase space we associate two fundamental frequencies (ν 1 , ν 2 ) ∈ [0, 1) 2 and display them in a two-dimensional frequency space. This is numerically done using a Fourier-transform-based frequency analysis [50,86,87]. For a given initial condition (x 0 , y 0 , z 0 ) we use 4096 points of the orbit and compute the fundamental frequencies ν 1 and ν 2 using two complex signals cos y n + i sin y n and x n + iz n respectively. These correspond to the motion in y-direction and the projection on the (x, z)-plane, respectively.
To decide whether an orbit is regular or chaotic, two consecutive segments of 4096 points are compared, giving frequencies (ν 1 , ν 2 ) and (ν 1 ,ν 2 ). As chaos indicator we use For a regular orbit, the frequencies of the consecutive segments are very similar so that δ is very small. On the other hand, for a chaotic orbit, the frequencies of successive segments will be very different. Thus an orbit is accepted as regular if δ < δ reg . It is important to point out that this analysis is based on finite-time information. Hence, some orbits which are accepted as regular 2d tori may actually be chaotic orbits. However, this is a common limitation of any tool for chaos detection, see e.g. Ref. [88].
For the central regular region we use 10 8 random initial conditions with x ∈ [π, 2π], z ∈ [π/2, 3π/2] and y ∈ [0, 2π[ and use δ reg = 10 −9 . This gives a total of 2.3 × 10 6 regular tori for C = 0.06. These frequency pairs are shown in the two-dimensional frequency space in Fig. 7. They are essentially arranged along a one-dimensional curve, emanating from the point of the frequencies of the central 1d torus at (ν c 1 , ν c 2 ) ≈ (0.4495, 0.325). This illustrates that the regular orbits form a one-parameter family in the volume-preserving case [62,89]. Note that in contrast in 4d symplectic maps regular orbits occur as two-parameter families, leading to two-dimensional regions in frequency space [90].
Of particular importance for the dynamics are so-called resonances for which the frequencies fulfill where m 1 , m 2 , n ∈ Z without a common divisor and at least one of m 1 or m 2 is not zero. The order of a resonance is given by |m 1 | + |m 2 |. These resonance lines, denoted by m 1 : m 2 : n, form a dense resonance web in frequency space. In Fig. 7, the curve of regular 2d tori in frequency is interrupted by several gaps. These gaps are related to resonance lines for which the corresponding original regular 2d tori are destroyed. The number of independent resonance conditions for a given frequency pair (ν 1 , ν 2 ) constitute the rank of the resonance and correspond to different types of dynamics: • If the frequency pair fulfills no resonance conditions, it is of rank-0. The motion on the corresponding 2d torus is quasi-periodic filling it densely. This corresponds to the dynamics on the regular cylinders illustrated in Fig. 2.
• If only one resonance condition is satisfied, the resonance is either (a) uncoupled, i.e., m 1 : 0 : n or 0 : m 2 : n, or (b) coupled, i.e., m 1 : m 2 : n with both m 1 and m 2 nonzero. Such rank-1 resonances correspond to quasi-periodic motion on a 1d invariant set which consists either of one component in the case of coupled resonances or of m 1 (or m 2 ) dynamically connected components in the case of uncoupled resonances. We would like to point out that this is the most important type of resonance for stickiness considered in the next section.
In Fig. 7 orbits belonging to the 2 : 0 : 1 resonance, the 0 : 4 : 1 resonance, and the 1 : 5 : 2 resonance, respectively, are shown. Note that the 2 : 0 : 1 resonance provides an example of regular orbits which do not fulfill the slice condition (4), and are therefore not visible in Fig. 4.
• For rank-2 resonances, two independent resonance conditions are fulfilled simultaneously. While there are a few examples of such double resonances for the ABC map, they appear to be of no relevance for the stickiness of orbits considered in the next section.
Due to the resonances, which form a dense set of lines in frequency space, the one-parameter family of 2d tori has infinitely many holes. Typically, resonances of low order lead to large gaps, such as the 2 : 0 : 1 resonance. To the right and below the frequencies (ν bc 1 , ν bc 2 ) of the boundary cylinder, all regular tori have to be resonant, as any non-resonant 2d torus would form a larger boundary cylinder. Numerically this is not easy to confirm, in particular for large orders. Moreover, as mentioned above, there could also be orbits for which the frequency criterion (5) is fulfilled, but which escape to the chaotic region for very long times.

III. STICKINESS AND POWER-LAW TRAPPING
To characterize transport and stickiness in a mixed phase space, a powerful approach is based on the Poincaré recurrence theorem.
It states that for a measure-preserving map with invariant probability measure µ, almost all orbits started in a region Λ in phase space will return to that region at a later time [91]. It is numerically convenient to consider the cumulative Poincaré recurrence statistics defined as where N (0) is the number of orbits initially started in Λ and N (t) is the number of orbits which have not yet returned to Λ until time t > 0. By definition, P (0) = 1 and P (t) decreases monotonically with time.
For the determination of the Poincaré recurrence statistics, we choose N (0) = 10 12 initial conditions randomly from a uniform distribution in Λ. Figure 8 shows the Poincaré recurrence statistics P (t) for the system at C = 0.06. An overall power-law decay is found, with several significant oscillations superimposed. We will explain their origin in Sec. III D. Note that the same type of behavior is also found for other values of C ∈ [0, 1] (not shown). The origin of the slow decay of P (t) are chaotic orbits which spend long times at either of the two regular regions. As both regular regions are related by symmetry, the contributions to P (t) should statistically be the same. To confirm this, we classify all long-trapped orbits to belong to either of the regions using the average value of their x coordinate. Both types of orbits give identical contributions to P (t) for t > 10 4 , while initial differences are due the non-symmetric location of Λ. Therefore, in the following analysis of long-trapped orbits, we restrict to those orbits which belong to the central region.

B. Long-trapped orbits
Long-trapped orbits are found to stay close to the surface of the boundary cylinder which constitutes an impenetrable barrier. As representative example we consider in the following the longest trapped orbit found in our numerical studies with t rec = 3 × 10 8 . Figure 9 shows the last 50 000 iterates of this orbit in the 3d phase space with time encoded in color. It fills a cylindrical surface outside of the boundary cylinder before moving away to the chaotic sea, as seen by the red scattered dots. The approach to the boundary cylinder is more clearly seen in the 2d slice by plotting all points of the longtrapped orbit fulfilling the slice condition (4), see Fig. 10. The magnification in the inset reveals that a small region is left out, which is due to the 5 : 15 : 7 resonance. In the 3d phase space the corresponding elliptic 1d torus completes 15 rotations on the x − z plane around the central 1d torus while it goes around 5 times in the y-direction. This therefore leads to 15 small resonance islands in the 2d slice.
To determine the boundary cylinder, shown as black curve in the 2d slice in Fig. 4 and in Fig. 10, we consider several initial conditions along a line parallel to the zaxis starting from the central 1d torus, see Fig. 10. For each initial condition the frequencies (ν 1 , ν 2 ) for the first 4096 points and (ν 1 ,ν 2 ) for 4096 points after 10 9 iterates are computed. The last initial condition for which the frequency criterion (5) is fulfilled provides a good approximation to the boundary cylinder, for which we get (ν bc 1 , ν bc 2 ) ≈ (0.5399, 0.2868). Note that this numerical approach of course cannot guarantee that this really is a 2d torus as for very large number of iterations it could eventually escape towards the chaotic region.
For a long-trapped orbit, such as shown in Fig. 9 and Fig. 10, the motion is close to the boundary cylinder. Therefore one can associate finite-time frequencies by computing a sequence of frequency pairs (ν k 1 , ν k 2 ), each  determined from 4096 consecutive points of the orbit. As the regular tori form a 1d line in frequency space, see Fig. 7, the approach to the boundary cylinder can be quantified by either ν 1 or ν 2 . Here we use ν 1 which corresponds to the motion in y-direction. Figure 11 shows for the long-trapped orbit the frequencies line 5 : 15 : 7 is shown. Crossing the 20 : −3 : 10 resonance, the orbit rapidly moves towards the frequency of the boundary cylinder and stays around the resonance line for most of the time. The same orbit is shown in frequency space in Fig. 12 with points colored according to time. While the initial (blue) and final (orange) part extend over some more scattered range in both ν 1 and ν 2 , most of the time is spend close to (ν bc 1 , ν bc 2 ). In addition several gaps are visible which are related to resonance lines, as we will discuss in more detail in the next section.
Note that to associate a frequency ν 1 with a resonance line we use the fact that the regular 2d tori form an essentially one-dimensional curve in frequency space, as seen in Fig. 7. Thus this line intersects a given resonance line m 1 : m 2 : n only in one point which provides a unique ν 1 , as displayed in Fig. 11. Numerically we fit the data for the regular 2d tori by a quadratic polynomial which is then used to compute the intersection points.

C. Approach to the boundary cylinder
Based on the set of long-trapped orbits determined using the Poincaré recurrence statistics we now analyze their statistical properties. In particular we identify the most relevant resonance lines and then use these to split the contributions to P (t) to explain the overall power-law behavior and the superimposed oscillations.

Identification of resonance lines
To identify the most relevant resonance lines, a set of all resonance lines up to order 25 is generated in the rectangle ν 1 ∈ [0.51, 0.59], ν 2 ∈ [0.26, 0.30] in frequency space. For one given long-trapped orbit the sequence of frequency pairs (ν k 1 , ν k 2 ) is computed and for each frequency pair the distance to a resonance line m 1 : m 2 : n is determined by The resonance with the smallest distance less than a threshold (set here as 10 −9 ) is associated with the k-th orbit segment. The resulting six most dominant resonance lines with which most of the frequency pairs of all orbits are associated are -9 : 9 : 8, 1 : 20 : 6, 5 : −3 : 2, −1 : 2 : 0, −1 : 9 : 2, and 20 : −3 : −10. Interestingly, all of them are coupled resonances. These lines are shown in frequency space in Fig. 12. When chaotic orbits approach the boundary cylinder, these most significant resonances are successively accessed. The approach is characterized by the distance to the boundary cylinder,   is already trapped with smaller ∆ν. In addition, a higher order resonance 10 : 37 : 16 appears to be relevant only for orbits t rec ∈ [10 7 , 5 × 10 7 ] in Fig. 13(g). In fact, the longest-trapped orbit shown in Fig. 11 is initially trapped around this resonance for about 10 7 iterations before it approaches the 5 : 15 : 7 resonances. Figure 14 shows the frequency ν 1 as function of time for two examples illustrating the relevance and geometry of coupled resonances. In Fig. 14(a) an orbit with t rec = 2.2 × 10 7 is shown, which initially approaches ν bc 1 and fluctuates around the frequency of the 20 : −3 : 10 resonance. Finally, the orbits sticks on the coupled −1 : 2 : 0 resonance for an even longer time. The corresponding dynamics in the 3d phase space is illustrated in the inset and consists of two separate, but dynamically connected parts. In Fig. 14(b) an orbit with return time t rec = 7.3 × 10 7 is shown. Here 5 : −3 : 2 is the dominant resonance. The corresponding dynamics in the 3d phase space shown in the inset spirals five times around in the y-direction and three times around the central 1d torus. For both orbits ν 1 shows on this scale essentially no fluctuations around the frequency of the corresponding resonance. A closer inspection using the representation in the 2d slice shows no indication that higher order resonances deeper in the class hierarchy around the corresponding resonance islands are accessed. Still we expect that this is possible for even longer trapped orbits.

D. Splitting of P (t)
The plot in Fig. 13(h) shows the complete distribution of ∆ν for all long-trapped orbits for t rec ≥ chambers as displayed in Fig. 15. The aim is to quantify the contributions to the Poincaré recurrence statistics P (t) from each of these chambers. For this we associate, based on ∆ν, with each long-trapped orbit the most relevant chamber, i.e. the one for which it spends the longest time. By this the total P (t) is split into contributions P i (t), determined from those orbits associated with the i-th chamber. Figure 16 shows these contributions together with the full P (t). Several of the P i (t) show an initially exponential decay. In contrast, for P 2 (t) an overall power-law is found for t 10 5 . This suggests that there is some generalized hierarchy (e.g. involving the equivalent of the class hierarchy in 2d maps [41]), which is accessed by orbits in this chamber. Moreover, chamber 2 also contributes significantly in the initial decay up to t 7 × 10 4 of the full P (t) shown in red. From t 7 × 10 4 on the most relevant contribution comes from chamber 3 for which P 3 (t) shows good agreement with an exponential decay, as manifested in a semilogarithmic plot giving a straight line over several orders of magnitude (not shown). Starting from approximately t = 3 × 10 6 one sees that P 8 (t) gives the most important contribution, again well-described by an exponential. In particular, the sum of the four most important P i (t) provides a good approximation of the full P (t), see Fig. 17. This therefore on the one hand explains the overall power-law behavior of P (t) as sum of exponentially decaying contributions, in the spirit of the explanation of power-law trapping in 2d systems [16]. On the other hand, a power-law requires a particular scaling of the individual contributions and is only expected to hold asymptotically for large t. The observed superimposed oscillations can be traced back to the individual contributions P i (t). In particular P 3 (t) is responsible for  the prominent hump near t = 10 6 . In general it has been shown that uncorrelated fluctuations in the transition rates between regions separated by partial barriers lead to slowly decaying oscillations in P (t) [31].

E. Partial barriers
Overall, the mechanism of trapping in the 3d ABC map appears to be similar to the well-known trapping in 2d maps: Long-trapped orbits approach a boundarycylinder and longer trapping times correspond to a closer approach, which is similar to the level hierarchy [92]. Also trapping within resonances is found, which is similar to the class hierarchy [41]. In the 2d case these hierarchies can be arranged in a binary tree [40,93], which under certain statistical assumptions leads to a universal asymptotic decay [30,34]. However, an important difference to the 2d case is that for the 3d ABC map coupled resonances are possible and appear to play the key role for long-trapped orbits. The numerics suggests that there must be some kind of partial barriers which define the borders between the chambers containing the identified resonances. Geometrically such partial barriers should correspond to broken cylindrical 2d tori. Here 2d tori which are furthest away from resonances should persist longest with increasing perturbation from the uncoupled case. The break-up of invariant tori in volume-preserving maps has been studied recently in [65,68] showing that tori seem to stretch near break-up which might lead to analogues of gaps in cantori for 2d maps. Computing a sequence of 1d tori approximating such partial barriers and determining the transport through them for volumepreserving maps such as the ABC map, is an interesting non-trivial task left for future studies.

IV. SUMMARY AND OUTLOOK
In this paper, we study stickiness in volume-preserving maps using the example of the ABC map. The 3d phase space of the map is mixed with two disjoint regular regions embedded in a chaotic sea. The map reduces to a (2+1)d system when C = 0. In this case the motion in the y-direction decouples from the 2d area-preserving x − z subsystem, which has two elliptic fixed points. For the uncoupled case regular regions of the 3d map are in the form of concentric cylinders organized around 1d elliptic tori corresponding to the fixed points of the 2d subsystem. With increasing C > 0 these cylinders are deformed. As long as these deformed cylinders acquire no holes, they form absolute barriers to the motion. If a cylinder is destroyed, this is expected to lead to a partial barrier to transport for chaotic orbits. For a finite C, destruction of 2d tori around 1d tori may occur inside the regular regions in addition to the interface of the chaotic and regular region outside. Thus chaotic layers appear inside the regular regions which are not yet connected to the chaotic sea, as is nicely seen in the sequence of plots of the finite-time Lyapunov exponents in Fig. 5. Representing the 2d regular tori in frequency space shows that they lead to an essentially one-dimensional curve which is attached to the frequencies of the central elliptic 1d torus. Resonances lead to gaps within the family of 2d tori. An outermost regular torus is identified as boundary cylinder which provides the separation to the chaotic sea.
Using the cumulative Poincaré recurrence statistics P (t) allows to investigate the stickiness of chaotic orbits. For the ABC map we find an overall power-law decay with significant oscillations superimposed. Longtrapped orbits approach the boundary cylinder and their segment-wise frequency-time representation reveals that they spend long times near six most relevant resonances. Using the distribution of all frequencies of all longtrapped orbits shows that they can be classified into nine chambers. Each long-trapped orbit is classified according to the chamber in which it spends the largest amount of time. From this we obtain a splitting of the cumulative Poincaré recurrence statistics P (t) into different contributions P i (t). These show good agreement with exponential decays and the sum of only four of them gives an excellent approximation to the full P (t). This explains both the overall approximate power-law decay, as a superposition of individual exponential decays, and also the superimposed oscillations due to the individual P i (t). This also suggests the existence of some kind of partial barriers separating the different chambers. Most likely, these are broken 2d tori, similar to cantori in 2d maps.
An important task for the future is to explicitly determine the partial barriers, e.g. by sequences of approximating families of 1d tori, and to identify the turnstile and to determine the corresponding flux. This should eventually allow for setting up a quantitative statistical prediction of the observed power-law trapping. Moreover, a fundamental question concerns whether there is a universal value of the exponent γ of the power-law decay, similar to that in area-preserving maps. In addition, power-law trapping indicates anomalous diffusion, so that a detailed study of the diffusive properties in the ABC map would be very interesting.