Areapreserving maps models of gyroaveraged
chaotic transport
Abstract
Discrete maps have been extensively used to model 2dimensional chaotic transport in plasmas and fluids. Here we focus on areapreserving maps describing finite Larmor radius (FLR) effects on chaotic transport in magnetized plasmas with zonal flows perturbed by electrostatic drift waves. FLR effects are included by gyroaveraging the Hamiltonians of the maps which, depending on the zonal flow profile, can have monotonic or nonmonotonic frequencies. In the limit of zero Larmor radius, the monotonic frequency map reduces to the standard ChirikovTaylor map, and, in the case of nonmonotonic frequency, the map reduces to the standard nontwist map. We show that in both cases FLR leads to chaos suppression, changes in the stability of fixed points, and robustness of transport barriers. FLR effects are also responsible for changes in the phase space topology and zonal flow bifurcations. Dynamical systems methods based on recurrence time statistics are used to quantify the dependence on the Larmor radius of the threshold for the destruction of transport barriers.
I Introduction
Important research efforts in controlled nuclear fusion are focused on the magnetic confinement of hot plasmas. In order to improve the confinement conditions, a better understanding of the particle transport is needed, especially in the case of particle drift motion. A standard approach to this problem is based on the approximation of charged particle’s guiding center’s motion, e.g. Refs. klevadrake1984 (); Pettini (); ibere_horton_2008 (); delcastillo2000 (). However, in the case of fast particles, e.g. alpha particles in burning plasmas, or in the presence of inhomogeneous fields on the scale of the Larmor radius, it is necessary to consider finite Larmor radius (FLR) effects to correctly estimate the transport of particles Manfredi96 (); delcastillomartinell2013 ().
Previous studies on the role of the Larmor radius on transport include Refs. Annibaldi02 (); Manfredi96 (); Manfredi97 () where particle transport in numerical simulations of electrostatic turbulence was analyzed and FLR effects were shown to inhibit transport. The problem of nondiffusive chaotic transport and fractional diffusion in the presence of FLR effects was addressed in Ref. gustafson (). In Refs. delcastillomartinell2012 (); delcastillomartinell2013 () dynamical systems methods were used to investigate FLR effects on chaos suppression and bifurcations in the phase space topology.
In this paper we study FLR effects in the context of simplified areapreserving map models of motion. The maps are constructed following the standard Hamiltonian framework for electrostatic drift motion horton1985 (), and FLR corrections are included by gyroaveraging the electrostatic potential Lee87 (). The potential has an equilibrium depending on the radial coordinate and a timedependent perturbation consisting of a linear superposition of drift waves. Depending on the radial dependence of the equilibrium the maps’ frequencies can have monotonic or nonmonotonic profiles. Our main focus is the study of FLR effects on chaos suppression, stability of fixed points, nontwist phase space topologies, and transport barriers. Although drift might play a role in the presence of FLR, here we will restrict attention to drifts.
In dynamical systems theory, the term “nontwist” is used to designate Hamiltonian systems that violate the twist or nondegeneracy condition. A paradigmatic example is a perturbed system in which the frequency of the unperturbed part of the Hamiltonian is a nonmonotonic function of the action. Nontwist Hamiltonian systems are found in many physical systems including transport in magnetized plasmas delcastillo2000 (); ibere_horton_2008 (), magnetic fields with reverse shear in toroidal plasma devices Oda (); Balescu98 (); delcastillo1992 (); Firpo2009 (), and transport by traveling waves in shear flows delcastillo1993 (); BeronVera10 (); Budyansky09 () among others. The presence of nontwist transport barriers is among the most important properties of nontwist Hamiltonian systems. By nontwist transport barriers we mean a robust region of invariant circles, also called KolmogorovArnoldMoser (KAM) curves, that are very resilient to breakup, i.e., they can survive even when the phase space is almost completely chaotic. This property, originally discovered in Refs. delcastillo1993 (); delcastillonegreteetal1996 (), has also been referred to as strong KAM stability Rypina () in the context of oneandahalf degreesoffreedom Hamiltonian systems, and has been the subject of several studies (see, for example, Portela07 (); Szezech09 (); Caldas (); delaLlave2000 () and references therein). Throughout this paper use the term “KAM curve” as a synonymous of “invariant circle” or “transport barrier”.
In this paper the maps correspond to Hamiltonians whose unperturbed frequencies can be monotonic or nonmonotonic functions of the action variable. The latter case is characterized by the presence of nontwist transport barriers that can be destroyed or restored by changing the value of the Larmor radius. This property is directly related to the FLR suppression of chaotic transport studied in Refs. delcastillomartinell2012 (); delcastillomartinell2013 (). Here we show that chaos suppression occurs when the Larmor radius is close to specific values for which invariant circles become very resilient to breakup. In the case of nontwist maps, the invariant circles forming the nontwist transport barrier are the easiest to restore and the hardest to break. One of the main goals of the present work is the numerical computation of the critical parameters for the breakup of nontwist barriers. To compute the breakup diagrams we follow the technique in Ref. Altmann () which is based on the application of Slater’s theorem Slater () to the recurrence properties of orbits. As shown in Refs. Altmann (); Zou (), the recurrence properties of an orbit provide an efficient and accurate method to differentiate chaotic and nonchaotic motion. A recent example of the application of this technique to the standard nontwist map was presented in Ref. Abud ().
The rest of the paper is organized as follows. The next section introduces a general area preserving map of gyroaveraged transport. This model is the starting point for the specific cases studied in the subsequent sections. The first case, for which the frequency has a monotonic profile, is presented in Sec. III, where FLR effects on the transition to global chaos are analyzed. Section IV describes a case with a nonmonotonic profile and discusses FLR effects on the stability of fixed points, phase space topology, and breakup diagrams. A second nontwist map is presented in Sec. V where FLR effects on zonal flow bifurcations are analyzed. Section VI presents the conclusions and a brief summary of the main results.
Ii Transport Model
The drift velocity of the guiding center is given by Nicholson ()
(1) 
Using as radial coordinate, and as poloidal coordinate, the equations of the drift motion, , can be written as the Hamiltonian system
(2) 
where
(3) 
is the electrostatic potential, and denotes the magnitude of the toroidal magnetic field, . As discussed in Ref. delcastillomartinell2012 (), finite Larmor radius (FLR) effects can be incorporated substituting the electrostatic potential by its average over a circle around the guide center
(4) 
where is the Larmor radius. Equation (4) corresponds to the wellknown gyroaveraging operation Lee87 (). The gyroaveraged Hamiltonian can then be defined as
(5) 
and the gyroaveraged equations of motion (2) can be written as
(6) 
Following Ref. horton1998 (), we assume an electrostatic potential of the form
(7) 
where is the equilibrium potential, , the amplitude of the drift waves, is the wave number, and is the fundamental frequency. Applying the gyroaverage operation in Eq. (4) to Eq. (7), and substituting the result in Eq. (3), we obtain the Hamiltonian
(8) 
where is the zeroorder Bessel function, and the integrable Hamiltonian is defined as
(9) 
Using the Fourier series representation of the Dirac delta function, Eq. (8) can be rewritten as
(10) 
Let and , with , , and . Integrating Eq. (6) in the interval leads to the gyroaveraged drift wave map
(11)  
(12) 
where , and corresponds to the frequency associated to the integrable Hamiltonian in Eq. (9)
(13) 
For , , where is the radial component of the electric field. Depending on the equilibrium potential , which determines and , different areapreserving maps can be constructed from Eqs. (11) and (12). In the next sections, we discuss three different cases.
Iii Gyroaveraged standard map
As a first simple example we assume a monotonic linear frequency profile to construct the gyroaveraged standard map (GSM). This map is a modified version of the standard map, also known as the ChirikovTaylor map chirikov79 (); taylor (). To define the frequency in Eq. (12), we use the equilibrium potential
(14) 
where is the wave number in Eq. (7), and is a free parameter. Applying the gyroaverage operation to Eq. (14) yields
(15) 
Substituting Eq. (15) in Eq. (13), we get the frequency
(16) 
which does not depend on the Larmor radius, and has a monotonic profile. Introducing the nondimensional variables
(17) 
and the constant
(18) 
we get, from Eq. (11) and Eq. (12), the gyroaveraged standard map (GSM)
(19)  
(20) 
where
(21) 
is the effective perturbation parameter and is the perturbation parameter. For , , and the GSM reduces to the standard map.
The phase space of the GSM, as in the case of any areapreserving map, consists of periodic, quasiperiodic, and chaotic orbits. Quasiperiodic orbits cover densely invariant curves. The invariant curves around the elliptic fixed points form island chains, and the invariant curves that wind around the entire domain of the angle variable form invariant circles Meiss92 (). The presence of KAM curves is of special interest because they preclude transport in the direction of the radial coordinate , keeping chaotic orbits confined to specific regions of phase space.
According to Greene’s residue method Greene79 (), the transition to global chaos in the standard map occurs when the absolute value of the perturbation parameter is equal to, or greater than, the critical value . That is, for , all KAM curves are broken, and chaotic orbits can spread over all phase space (except in regions occupied by isolated islands). Larmor radius effects on the transition to global chaos in the GSM can be analyzed by defining the critical line dividing the parameter space in two regions: one for which the phase space contains at least one KAM curve and another for which transition to global chaos has occurred. In the standard map, there is no dependence on , which means that the critical line is just a horizontal line defined at , as indicated in Fig. 1. However, for the GSM, the critical line is determined by the condition which implies
(22) 
As shown in Fig. 1, even for high values of the perturbation parameter () there are an infinite number of regions of Larmor radius values for which KAM curves can be restored and chaos is suppressed. The gyroaveraging operation “breaks” the critical line at the zeros of the zeroorder Bessel function. In particular, near the zeros, the critical perturbation goes to infinity and the transition to global chaos cannot occur. Writing , where is the th zero of and , and expanding in Taylor series, the divergence of the critical perturbation near the zeros of can be approximated as
(23) 
where denotes the Bessel function of first order. The first five positive zeros, indicated in Fig. 1, are approximately , , , , and .
The effect of chaos suppression is illustrated by Figs. 2 (a)(b). Figure 2(a) shows a GSM Poincaré section for and , which correspond to the point in Fig. 1. No KAM curves are observed in the Poincaré section as belongs to the region of global chaos. Keeping the same value of and changing the value of to (point in Fig. 1) eliminates all the chaotic orbits (see figure 2(b)).
Iv Gyroaveraged standard nontwist map
As a second example of the gyroaveraged drift wave map, we introduce the gyroaveraged standard nontwist map (GSNM), which corresponds to a nonmonotonic radial electric field. In this case the equilibrium potential is
(24) 
where and are dimensional constants. The gyroaverage of Eq. (24) gives
(25) 
and from Eq. (13) we get the nonmonotonic parabolic frequency profile,
(26) 
which has a maximum at . For , the zeros of (and also of the profile) are located at , where is the characteristic length of the frequency profile. Substituting Eq. (26) into Eq. (12), we obtain
(27)  
(28) 
Introducing the dimensionless variables:
(29) 
we get the gyroaveraged standard nontwist map (GSNM):
(30)  
(31) 
where:
(32)  
(33) 
are four dimensionless parameters. We refer to as the perturbation parameter, which is proportional to the amplitude of the drift waves. Like in the GSM case, we can also define an effective perturbation parameter
(34) 
As indicated in Eq. (33), , and correspond to the Larmor radius normalized using two different length scales: the characteristic length of the frequency profile and the wavelength respectively
iv.1 Fixed Points
According to Eqs. (30)(31), the fixed points of the GSNM, and , satisfy
(35)  
(36) 
where is an integer number. For each , , , and , there are four fixed points:
(37)  
(38) 
where
(39) 
Note that the coordinate of and does not depend on any parameters. If , the pair of points collide at , and collide at . Figure 3 shows the axis position of both pairs for and increasing . The collision occurs for . For higher values of , the fixed points don’t exist.
The stability of a periodic orbit of a map is determined by the residue
(40) 
where is the Jacobian matrix of , evaluated at the periodic orbit . If , the periodic orbit is elliptic (stable); if or , it is hyperbolic (unstable); and it is parabolic for or . Applying formula (40) to the GNSM fixed points in Eqs. (37)(38) we get
(41) 
The stability of the fixed points with can be analyzed using the parameter . As shown in Fig. 4, depending on the value of , there are three possible configurations. The symbol “x” denotes an hyperbolic point, and “o” an elliptic point. The points have their stability inverted when (configuration II), as indicated in Fig. 4(b). All the fixed points are hyperbolic in configuration III, Fig. 4(c), which occurs for or .
iv.2 Separatrix reconnection
The location and stability of the fixed points determine the different phase space topologies of the GSNM. These topologies, illustrated in Fig. 5, are characteristic of nontwist maps and are called heteroclinic, separatrix reconnection, and homoclinic delcastillo1993 (). Since the Larmor radius changes the stability of the fixed points, it is expected that it will also change the topology. Figures 5(a)(c) show the heteroclinictype, separatrix reconnection, and the homoclinictype topologies in Poincaré sections of the GSNM for different values of and , and fixed values of and .
To determine the condition for separatrix reconnection associated to the fixed points and with , we follow Refs. delcastillonegreteetal1996 (); delcastillo1993 () and approximate the GSNM in the vicinity of the fixed points by the Hamiltonian
(42) 
For (configuration I), separatrix reconnection occurs when . If (configuration II), reconnection is observed when . Combining these two conditions, we obtain the reconnection line
(43) 
with slope:
(44) 
The reconnection line divides the  parameter space in two regions: one with heteroclinictype topology () and another with homoclinictype topology (). The slope of the reconnection line is defined by the angle . Figure 6 shows isolines of the angle . Keeping and fixed, the topology of the phase space remains unchanged if the parameters and vary over an isoline. As an example, for and such that (which corresponds to the reconnection condition in the standard nontwist map delcastillo1993 ()) there is reconnection for any values of and over the red isoline in Fig. 6. The red isoline crosses the point where the GNSM reduces to the standard nontwist map. For high (red) values of , the slope of the reconnection line approaches , and the phase space is characterized by the homoclinictype topology for a fixed and almost every value of the parameter . For low (blue) values of , the slope of the reconnection line tends to , and the phase space is characterized by the heteroclinictype topology for a fixed and almost every value of the parameter .
iv.3 Nontwist transport barrier
The transition to global chaos in the GSNM corresponds to the destruction of the nontwist transport barrier (NTB) which is a robust nonchaotic region of KAM curves dividing the phase space. The NTB is robust in the sense that their KAM curves are typically more resilient to perturbation than other KAM curves in the system.
Consider the degreesoffreedom Hamiltonian system
(45) 
where and are the actionangle variables of the integrable system, , and the parameter controls the strength of the timeperiodic perturbation . In this case, the twist (nondegeneracy) condition is
(46) 
for all . If there is at least one critical point for which , the nondegeneracy condition is violated the Hamiltonian in Eq. (45) is nontwist. One of the main properties of these systems is the presence of robust NTBs located at the critical points (maxima or minima) of the frequency profile even in the case when most of the phase space is chaotic delcastillo1993 (). As discussed before, the GSNM can be obtained by integrating the equations of motion of the degrees of freedom Hamiltonian in Eq. (8). This is a nontwist Hamiltonian because the frequency in Eq. (26) has a critical point at . Figure 7 shows a Poincaré section of the GSNM with , , and . As expected, the phase space is characterized predominantly by chaotic orbits (green and red regions), but it has a robust NTB consisting of “belt” of KAM curves (black curves).
iv.4 Breakup Diagrams
Here we study the transition to global chaos in the GSNM, i.e., the destruction of the robust KAM curves that constitute the NTB. It is important to remark that the transition to global chaos in nontwist systems is still an open problem. However, a possible approach to this problem consists of estimating the parameter values, or the critical thresholds, that determine the breakup of the shearless KAM curve. In the absence of a perturbation, the shearless KAM curve is located where the twist condition is violated, and the parameter values for the destruction of the NTB can be estimated from the parameter values that determine the breakup of the shearless curve.
A standard way to determine the parameter values for the destruction of the shearless curve is based on the indicator point (IP) method proposed in Ref.shinohara () and used is several works including Refs. Wurm (); Budyansky09 (); delcastillomartinell2013 (). This method consists on finding the parameter values for which the iterations of the indicator point, which defines the indicator point orbit (IPO), are chaotic. If the IPO is quasiperiodic, it traces a shearless KAM curve. When the IPO is chaotic, the shearless curve is broken. In nontwist maps with special symmetries, like the ones studied here, the indicator points can be easily found by computing the fixed points of the reversing symmetry group transformations. Further details can be found in Ref. Petrisor (). Following a procedure similar to the one in Ref. Wurm (), we obtained the following IPs for the gyroaveraged standard nontwist map
(47) 
(48) 
Once the IPs are found, the next step is to obtain the critical thresholds using an appropriate criterion to determine if the IPO is chaotic. Our criterion is based on the technique proposed in Altmann (), which considers the recurrence properties of the IPO in conjunction with Slater’s Theorem Slater (). Recurrence properties and their relation to the Slater’s Theorem were also discussed in Zou (). In Ref. Abud (), the technique described in Ref. Altmann () was used to study the breakup of the shearless curve in the nontwist standard map. The Slater’s Theorem states that, for any quasiperiodic motion on the circle, there are at most three different recurrences, or return times, in any connected interval. Although this theorem was originally formulated for circle maps, it can also be applied to twodimensional systems, once the dynamics on a KAM curve is mapped to a quasiperiodic rotation on a circle Altmann (). The recurrence time is defined as the number of iterations that an orbit takes to return to a neighborhood of a point. Our procedure to characterize the dynamics of an indicator point orbit of the GSNM is the following: with any one of the indicator points as the initial condition, we compute the orbit for iterations and choose the point of the orbit which has the maximum number of different recurrence times. If there are more than three different recurrence times, we conclude that the orbit is chaotic and the corresponding KAM curve is destroyed. More precisely, the procedure consists of the following steps:

Compute the orbit with initial condition , where is the number of iterations.

Construct the recurrence matrix ^{1}^{1}1Here we borrow the definition presented in Ref. Zou () of the binary matrix . Its graphical representation, called recurrence plot, can be used to analyze the recurrence properties of any dynamical system.:
(49) where , , is the Heaviside function, and is a parameter defining the size of the neighborhood. If the distance between and , given by the norm , is less than , ; otherwise, .

Define the recurrence time as for and . That is, the recurrence time is if the orbit crosses the neighborhood of at the th iteration. For each point , compute the set of different recurrences .

Determine the maximum number of different recurrence times:
(50) where is the number of recurrence times belonging to , and use the Slater’s Theorem to conclude that the orbit is chaotic if .
For , the characterization of the dynamics using the Slater’s Theorem is inconclusive because the orbit might be periodic, quasiperiodic or even chaotic. In this case, the number of iterations, , can be increased to explore if increases beyond . If the condition still persists and no change in is observed, the orbit is deemed periodic if and quasiperiodic if . The case requires a more careful analysis with higher values . Increasing , stabilizes in , indicating quasiperiodic dynamics, or assumes values greater than , which is the case for chaotic dynamics. Although rigorous criteria for the case are lacking, the numerical results in Refs. Altmann (); Zou () support the use of recurrences in conjunction with Slater’s theorem as a computationally efficient diagnostic to identify chaotic orbits.
We now apply this diagnostic to compute breakup diagrams for the GSNM showing regions in parameter space where the shearless curve is broken. The diagrams were constructed by computing the recurrences of IPOs for , and two different number of iterations, and . Since the GSNM has four free parameters only two parameters were varied at a time while the others were kept constant. Our first example is presented in Fig. 8 where we fixed , (which corresponds to the limit ), and plotted the values of on a grid in the space. The main difference between and is observed in points with , but even for a relatively small number of iterations it is possible to identify the breakup of the shearless curve in domains where . These domains have points plotted with magenta, red, and yellow colors. Green points () are concentrated in domains with low effective perturbation , which means small or close to a zero of the zeroorder Bessel function. For parameter values outside the green region (e.g., points and in Fig. 8) the IPO is chaotic, and, as shown in Fig. 9, the shearless curve is destroyed. For values near a zero of (e.g., point of Fig. 8), the shearless KAM curve is restored as shown by the red orbit in Fig. 9.
Breakup diagrams with varying and are shown in Fig. 10. The number of iterations used to construct the first diagram was . In the second diagram, the number of iterations was increased to . As before, regions where , which are indicators of nonchaotic dynamics, are concentrated near the zeros of . A high occurrence of points with is also observed in regions with low . When the number of iterations is increased from to , it s observed that the main geometric features remain the same, although several points with and transition to . Poincaré sections for parameters corresponding to points , , and in Fig. 10 are showed in Fig. 11ac. As expected, the shearless KAM curve is broken at points and , and restored at point .
Additional examples of breakup diagrams are presented in Figs. 12 (a)(d). The parameters are the same as those in Figs. 8(b) and 10(b), but the value of the fixed parameter is higher. Comparing Figs. 12 (a)(b) to 8(b), it is observed that increasing leads to an increase in the number of points with and to the robustness of the shearless KAM curve. As shown in Figs. 12 (c)(d), the distribution of points with changes when increases, resulting in different critical thresholds but with no significant suppression of green points. Red and yellow points disappear in certain regions and reappear in others, as can be seen by comparing Figs. 10(b), 12(c) and 12(d).
Summarizing, in this section we used a computational technique based on the recurrence properties of the IPO to estimate the critical parameter values for the breakup the shearless KAM curve. In particular, we computed breakup diagrams to understand the role of finite Larmor radius effects on the destruction and formation of shearless KAM curves. It was observed that when is close to a zero of , the shearless curve becomes more robust to perturbation in . This is because, near the zeros of , the GSNM’s Hamiltonian is effectively integrable. In particular, as approaches a zero of , the shearless curve, and the KAM curves that make up the nontwist transport barrier, are restored. The robustness of the shearless curve and the critical thresholds are also modified by increasing .
V Gyroaveraged quartic nontwist map
In this section, we propose another areapreserving map, the gyroaveraged quartic nontwist map (GQNM). A key property of the GQNM is that it exhibits a zonal flow bifurcation, similar to the one observed in the nontwist Hamiltonian system of Ref. delcastillomartinell2013 (). The Hamiltonian system in Ref. delcastillomartinell2013 () is a driftwave model of the transport with FLR effects, and the zonal flow bifurcation corresponds to a bifurcation of the maximum of the zonal flow velocity when the Larmor radius increases.
As in the previous maps, the GQNM is obtained from the gyroaveraged drift wave map in Eqs. (11) and (12). The equilibrium potential similar to the one proposed in Ref. delcastillomartinell2013 ()
(51) 
where and are dimensional constants. Applying the gyroaverage operation (4) to (51) gives
(52) 
and from Eq. (13) we get the nonlinear frequency
(53) 
where , and . The zonal flow bifurcation results from a bifurcation of the critical point of the frequency profile. To analyze it in a simple setting, consider the Taylor series approximation of (53) for small values of and ,
(54) 
Figure 13 shows the frequency profile in Eq. (54) for different values of . For , there is only one critical point, , which corresponds to a maximum at . Increasing leads to an increase of the profile’s “flatness” and eventually to a transition in which the maximum bifurcates into a minimum and two maxima. The bifurcation threshold can be determined from the condition,
(55) 
which for Eq. (54) gives
(56) 
Defining , , and , and using Eq. (54), we get from Eqs. (11) and (12), the GQNM
(57)  
(58) 
where , and the perturbation parameter is . As in the previous maps, an effective perturbation parameter can be defined as . Like the GSNM, the GQNM reduces to the standard nontwist map for . In what follows, we discuss FLR effects on the GQNM’s fixed points and nontwist transport barriers.
v.1 Fixed Points and nontwist transport barriers
The fixed points of the GQNM satisfy
(59)  
(60) 
For , , and , there are four fixed points independently of the value of ,
(61) 
If , there is an additional set of fixed points,
(62) 
where
(63) 
Figure 14 shows the coordinates of the fixed points with . and exist for all and are always located at . and only exist for and their positions in the axis are determined by the function , which satisfies .
The stability of the fixed points can be analyzed by evaluating the Green’s residue in Eq. (40)
(64)  
(65) 
where
(66) 
As shown in Fig. 15, depending on the value of , there are five possible stability configurations. The number of elliptic fixed points reduces when increases, and for all the fixed points are hyperbolic. In general, the first points to lose stability are the outer ones. Note that as approaches the zonal bifurcation threshold, diverges, and the range of the configurations III and IV increases.
An important consequence of the zonal flow bifurcation is the occurrence of nontwist barrier bifurcations. In particular, for , only one region with robust KAM curves is present (Fig. 16)(a), but for (Fig. 16)(b) two additional NTBs appear. The outer NTBs are separated from the central one by two regions of confined chaotic orbits.
v.2 Robustness of the central shearless curve
Our interest in studying the robustness of KAM curves as function of resides on the fact that this parameter is proportional to the amplitude of the drift waves. For very small Larmor radius (or in the absence of FLR corrections), high amplitude values can easily destroy all KAM curves. However, this is not always the case when FLR effects are taken into account. In particular, like in the GSNM (see Sec. IV.3), the robustness of the GQNM’s shearless curve is significantly increased for near a zero of .
Figures 17(a)(d) show breakup diagrams for the shearless KAM curve in the GQNM as function of and . In all cases, is small, because the values of are close to the zeros of . Under this condition, and when is close to the zonal flow bifurcation threshold in Eq. (56), a significant increasing of the robustness of the shearless curve is observed. For , higher values of are required for breakup as approaches , and for , the robustness is reduced as increases. That is, the robustness of the central NTB increases for close to the zeros of , and for close to .
It is interesting to note that, in the neighborhood of the shearless point, the degree of “flatness” of the frequency’s profile has a similar behavior. In particular, like the robustness of the central shearless curve, the degree of “flatness” (measured by ) is higher near the bifurcation threshold, but it becomes smaller elsewhere. These results provide support to the conclusion that the robustness of NTBs in nontwist maps depends on the flatness of the frequency profile at the maximum or minimum (critical points), and that, for small perturbation, the robustness increases with the degree of flatness.
The understanding of the role of the shear and the flatness of the frequency on the robustness of KAM curves is still an open problem. Among the first works addressing the role of shear is Ref. finn_1975 () where it was analytically and numerically shown that high shear reduces the strength of the perturbation required to break KAM surfaces. However this work limited attention to monotonic frequencies with nonvanishing shear (i.e., twist systems) and did not address the breakup of shearless KAM curves. The robustness of NTBs is also related to the strong KAM stability in nontwist maps Rypina (). According to the resonant overlap criterion chirikov79 (), tori located between resonances break when the resonances overlap. Thus, one expects tori to be more resilient when the resonances’ widths are small. As shown in Ref. Rypina (), in degreesoffreedom nontwist Hamiltonian systems, the width, , of secondorder degenerate resonances scales as
(67) 
where is the amplitude of the perturbation, and is such that
(68) 
which according to Eq. (46) corresponds to the point where the twist condition is violated. When the system is slightly perturbed, secondorder degenerate resonances appear in the neighborhood of the shearless curve, and their overlap leads to the break up of the shearless curve.
The central shearless curve of the GQNM is associated to the critical point where the frequency in Eq. (54) satisfies (68) for . Thus, in the Hamiltonian system from which the GQNM is derived, secondorder degenerate resonances might appear near the central shearless curve, and according to Eqs. (67) and (55), the width of these resonances scales as
(69) 
because vanishes for . Small values of result in small degenerate resonance widths and a more robust central shearless curve. As is directly related to the flatness of the frequency’s profile at , Eq. (67) is consistent with the results Figs. 17(a)(b) and supports the conjecture relating the robustness of shearless curves to the flatness of the frequency at the degenerate points. Similar arguments can be applied to interpret the results in Ref. Firpo2009 () where a study of the effect of the flatness of the profile on the robustness of NTBs was presented.
Vi Summary and Conclusions
In this paper we studied finite Larmor radius effects on discrete models of chaotic transport of test charged particles. The models are areapreserving Hamiltonian maps that assume an electrostatic potential consisting of the superposition of an equilibrium radial electrostatic field and a perturbation with a broad spectrum of drift waves. FLR effects were included by gyroaveraging the electrostatic potential and the maps were constructed from the integration of the resulting gyroaveraged Hamiltonian system.
We constructed and studied three maps: the gyroaveraged standard map (GSM), the gyroaveraged standard nontwist map (GSNM), and the gyroaveraged quartic nontwist map (GQNM). In the GSM, the frequency, i.e. the gyroaveraged zonal flow velocity, has a monotonic radial profile and the map satisfies the twist condition. In the GSNM and GQNM cases the frequency is nonmonotonic and as a result the maps are nontwist. Nonmonotonic frequencies have degenerate points (i.e., points where the frequency is maximum or minimum), and in the vicinity of these points robust nontwist transport barriers (NTBs) tend to form.
Typically, in area preserving maps, increasing the perturbation increases the chaos in the phase space. However, even for large perturbations, FLR effects can suppress chaos and restore broken KAM curves. The breakup diagrams provide the critical thresholds for the breakup of the shearless curve which is one of the robust KAM curves that constitute the NTB. To determine the breakup of the shearless curve we used an efficient and accurate method based on the maximum number of different recurrence times. We showed that a relatively low number of iterations is sufficient to determine the main features of the breakup diagrams. Although in general chaos is suppressed when the Larmor radius is close to zeros of the zeroorder Bessel function, the breakup diagrams exhibit a highly nontrivial fractallike dependence on the Larmor radius and the other parameters of the map.
We also studied the role of the Larmor radius on the topology of the phase space and on bifurcations of the zonal flow. In particular, we showed that the GSNM exhibits FLR dependent separatrix reconnection and obtained a formula determining the threshold for bifurcation from homoclinic to heteroclinic topologies as function of the Larmor radius. We also found that in the GQNM there is a critical value of the Larmor radius for which the maximum of the zonal flow bifurcates into a minimum and two maxima. These zonal flow bifurcations create additional shearless curves that play a highly nontrivial role on transport.
Motivated by the key role that periodic orbits play on transport we presented a detailed study of FLR effects on the location and stability of periodone orbits. In the GSNM there are four periodone orbits and, depending on the value of the Larmor radius, there are three possible stability configurations. In the GQNM the situation is more complex due to the existence of zonal flow bifurcations. In particular, there can be up to eight periodone fixed points, four of which only exist for Larmor radii greater than a critical threshold. In this case, depending on the value of the Larmor radius, there are five stability configurations. It was found that the elliptic fixed points created by varying the Larmor radius are strongly stable. The presence of elliptic points results in the formation of islands that can trap orbits in the phase space.
An important consequence of zonal flow bifurcations is the occurrence of up to three nontwist transport barriers. We focused on the central NTB and analyzed the robustness of the shearless curve using breakup diagrams. The results showed that the robustness of the central shearless curve significantly increases when the Larmor radius is close to the critical threshold for the zonal flow bifurcation. Based on this we argued that, in agreement with previous works, the “flatness” of the frequency profile at the critical point is directly associated with the robustness of the central shearless curve.
Vii Acknowledgments
This work was made possible through financial support from the Brazilian research agency FAPESP under grant 2013/004831. The authors would like to acknowledge Celso V. Abud for valuable discussions. DdCN was sponsored by the Office of Fusion Energy Sciences of the US Department of Energy at Oak Ridge National Laboratory, managed by UTBattelle, LLC, for the U.S. Department of Energy under contract DEAC0500OR22725. JDF gratefully acknowledges the hospitality of the Oak Ridge National Laboratory where part of work was conducted.
References
 (1) R. G. Kleva, and J. F. Drake, Phys. Fluids 27, 1686 (1984).
 (2) M. Pettini, A. Vulpiani, et al. Phys. Rev. A 38, 344 (1988)
 (3) F. A. Marcus, I. L. Caldas, Z. O Guimaraes Filho, P. J. Morrison, W. Horton, Y. K Kuznetsov and I. L. Nascimento, Phys. Plasmas 15, 112304 (2008).
 (4) D. delCastilloNegrete, Phys. Plasmas 7, 1702 (2000).
 (5) G. Manfredi and R. O. Dendy, Phys. Rev. Lett., 76, 4360, (1996).
 (6) J.J. Martinell and D. delCastilloNegrete, Phys. Plasmas 20, 022303 (2013).
 (7) S.V. Annibaldi, G. Manfredi, and R. O. Dendy, Phys. Plasmas 9, 791 (2002).
 (8) G. Manfredi and R. Dendy, Phys. Plasmas 4, 628 (1997).
 (9) K. Gustafson, D. delCastilloNegrete and W. Dorland, Phys. Plasmas 16, 102309 (2008).
 (10) D. delCastilloNegrete and J.J. Martinell, Commun. Nonlinear Sci. Numer. Simulat. 17, 2031 (2012).
 (11) W. Horton, Plasma Phys. Controlled Fusion 27, 9 (1985).
 (12) W. W. Lee, J. Comp. Phys 72, 243 (1987).
 (13) G. A. Oda and I. L. Caldas, Chaos, Solitons and Fractals, 5 15 (1995).
 (14) R. Balescu, Phys. Rev. E 58, 3 (1998).
 (15) D. delCastilloNegrete and P. J. Morrison, Bull. Am. Phys. Soc., Series II, 37, 1543 (1992)
 (16) L. Nasi and M.C. Firpo, Plasma Physics and Controlled Fusion, 51, 045006 (2009).
 (17) D. delCastilloNegrete and P. J. Morrison, Phys. Fluids A 5, 948 (1993)
 (18) F. J. BeronVera, M. J. Olascoaga, M. G. Brown MG, et al. Chaos 20, 017514 (2010).
 (19) M. V. Budyansky, M. Yu. Uleysky, and S. V. Prants, Phys. Rev. E, 79, 056215 (2009).
 (20) D. delCastilloNegrete, J. M. Greene, and P. J. Morrison, Physica D 91, 1 (1996).
 (21) I.I Rypina, M. G. Brown, et al. Phys. Rev. Lett. 98, 104102 (2007).
 (22) J. S. E. Portela, I. L. Caldas et al. International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, 17, 1589, (2007).
 (23) J. D. Szezech, I. L. Caldas, et al. Chaos: An Interdisciplinary Journal of Nonlinear Science, 19, 043108 (2009).
 (24) I.L. Caldas, R.L. Viana et al. Commun. Nonlinear Sci. Numer. Simulat. 17, 2021 (2011)
 (25) A. Delshams and R. delaLlave, SIAM J. Math. Anal., 31, 1235 (2000)
 (26) E. G. Altmann, G. Cristadoro, and D. Pazó, Phys. Rev. E, 73, 056201 (2006).
 (27) Y. Zou, D. Pazó, et al. Phys. Rev. E, 76, 016210 (2007).
 (28) C.V. Abud, PhD Thesis, University of São Paulo, (2013).
 (29) D.R Nicholson, “Introduction to Plasma Theory”. John Wiley and Sons, Inc. New York, 1983.
 (30) W. Horton, H.B. Park, J.M. Kwon, D. Strozzi, P. J. Morrison, and D.I. Choi, Phys. Plasmas 5, 11 (1998).
 (31) B. V. Chirikov, Phys. Rep., 52, 263, (1979).
 (32) J. B. Taylor, Culham Lab. Prog. Report CLMPR12, (1969).
 (33) J.D. Meiss, Rev. Mod. Phys. 64, 3 (1992).
 (34) J.M Greene, J. Math. Phys. 20, 1183 (1979)
 (35) S. Shinohara and Y. Aizawa, Prog. Theor. Phys. 97, 379 (1997).
 (36) E. Petrisor, Int. J. Bifurcation Chaos, 11, 497 (2001)
 (37) A. Wurm, A. Apte, et al. CHAOS 15, 2023108 (2005).
 (38) N. Slater, Proc. Cambridge Philos. Soc. 63, 1115 (1967).
 (39) J. M. Finn, Nuclear Fusion, 5, 845 (1975).