Area-preserving maps models of gyroaveraged {\bf E}\times{\bf B} chaotic transport

# Area-preserving maps models of gyroaveraged E×B chaotic transport

J. D. da Fonseca Physics Institute, University of São Paulo
São Paulo, SP, 5315-970, Brazil
D. del-Castillo-Negrete Oak Ridge National Laboratory
Oak Ridge, TN, 37831-8071, USA
I. L. Caldas Physics Institute, University of São Paulo
São Paulo, SP, 5315-970, Brazil
July 15, 2019
###### Abstract

Discrete maps have been extensively used to model 2-dimensional chaotic transport in plasmas and fluids. Here we focus on area-preserving 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 gyro-averaging the Hamiltonians of the maps which, depending on the zonal flow profile, can have monotonic or non-monotonic frequencies. In the limit of zero Larmor radius, the monotonic frequency map reduces to the standard Chirikov-Taylor map, and, in the case of non-monotonic 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. kleva-drake-1984 (); Pettini (); ibere_horton_2008 (); del-castillo-2000 (). 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 (); del-castillo-martinell2013 ().

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 non-diffusive chaotic transport and fractional diffusion in the presence of FLR effects was addressed in Ref. gustafson (). In Refs. del-castillo-martinell2012 (); del-castillo-martinell2013 () 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 area-preserving map models of motion. The maps are constructed following the standard Hamiltonian framework for electrostatic drift motion horton-1985 (), and FLR corrections are included by gyro-averaging the electrostatic potential Lee87 (). The potential has an equilibrium depending on the radial coordinate and a time-dependent perturbation consisting of a linear superposition of drift waves. Depending on the radial dependence of the equilibrium the maps’ frequencies can have monotonic or non-monotonic 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 non-degeneracy condition. A paradigmatic example is a perturbed system in which the frequency of the unperturbed part of the Hamiltonian is a non-monotonic function of the action. Nontwist Hamiltonian systems are found in many physical systems including transport in magnetized plasmas del-castillo-2000 (); ibere_horton_2008 (), magnetic fields with reverse shear in toroidal plasma devices Oda (); Balescu98 (); del-castillo-1992 (); Firpo2009 (), and transport by traveling waves in shear flows del-castillo-1993 (); Beron-Vera10 (); 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 Kolmogorov-Arnold-Moser (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. del-castillo-1993 (); del-castillo-negrete-etal-1996 (), has also been referred to as strong KAM stability Rypina () in the context of one-and-a-half degrees-of-freedom Hamiltonian systems, and has been the subject of several studies (see, for example, Portela07 (); Szezech09 (); Caldas (); de-la-Llave-2000 () 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 non-monotonic 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. del-castillo-martinell2012 (); del-castillo-martinell2013 (). 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 non-chaotic 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 gyro-averaged 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 non-monotonic 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 ()

 →vGC=→E×→BB2. (1)

Using as radial coordinate, and as poloidal coordinate, the equations of the drift motion, , can be written as the Hamiltonian system

 dxdt=−∂H(x,y,t)∂y,dydt=∂H(x,y,t)∂x, (2)

where

 H(x,y,t)=ϕ(x,y,t)B0, (3)

is the electrostatic potential, and denotes the magnitude of the toroidal magnetic field, . As discussed in Ref. del-castillo-martinell2012 (), finite Larmor radius (FLR) effects can be incorporated substituting the electrostatic potential by its average over a circle around the guide center

 ⟨ϕ(x,y,t)⟩φ=12π∫2π0ϕ(x+ρcosφ,y+ρsinφ,t)dφ, (4)

where is the Larmor radius. Equation (4) corresponds to the well-known gyro-averaging operation Lee87 (). The gyro-averaged Hamiltonian can then be defined as

 ⟨H(x,y,t)⟩φ=⟨ϕ(x,y,t)⟩φB0 (5)

and the gyro-averaged equations of motion (2) can be written as

 dydt=∂⟨H⟩φ∂x,dxdt=−∂⟨H⟩φ∂y. (6)

Following Ref. horton-1998 (), we assume an electrostatic potential of the form

 ϕ(x,y,t)=ϕ0(x)+A+∞∑m=−∞cos(ky−mω0t), (7)

where is the equilibrium potential, , the amplitude of the drift waves, is the wave number, and is the fundamental frequency. Applying the gyro-average operation in Eq. (4) to Eq. (7), and substituting the result in Eq. (3), we obtain the Hamiltonian

 ⟨H(x,y,t)⟩φ=⟨H0(x)⟩φ+AB0J0(kρ)+∞∑m=−∞cos(ky−mω0t), (8)

where is the zero-order Bessel function, and the integrable Hamiltonian is defined as

 ⟨H0(x)⟩φ=⟨ϕ0(x)⟩φB0 (9)

Using the Fourier series representation of the Dirac delta function, Eq. (8) can be rewritten as

 ⟨H(x,y,t)⟩φ=⟨H0(x)⟩φ+2πAB0J0(kρ)cos(ky)+∞∑m=−∞δ(ω0t−2πm). (10)

Let and , with , , and . Integrating Eq. (6) in the interval leads to the gyro-averaged drift wave map

 xn+1 (11) yn+1 =yn+2πω0Ω(xn+1) (12)

where , and corresponds to the frequency associated to the integrable Hamiltonian in Eq. (9)

 Ω(x)=d⟨H0(x)⟩φdx. (13)

For , , where is the radial component of the electric field. Depending on the equilibrium potential , which determines and , different area-preserving 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 gyro-averaged standard map (GSM). This map is a modified version of the standard map, also known as the Chirikov-Taylor map chirikov79 (); taylor (). To define the frequency in Eq. (12), we use the equilibrium potential

 ϕ0(x)=α(kx)22, (14)

where is the wave number in Eq. (7), and is a free parameter. Applying the gyro-average operation to Eq. (14) yields

 ⟨H0(x)⟩φ=αB0[(kx)22+^ρ24]. (15)

Substituting Eq. (15) in Eq. (13), we get the frequency

 Ω(x)=αk2B0x, (16)

which does not depend on the Larmor radius, and has a monotonic profile. Introducing the non-dimensional variables

 I=kγx,θ=ky, (17)

and the constant

 γ=2παk2ω0B0, (18)

we get, from Eq. (11) and Eq. (12), the gyro-averaged standard map (GSM)

 In+1= In+Kef(^ρ)sinθn (19) θn+1= θn+In+1,mod2π (20)

where

 Kef=KJ0(^ρ) (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 area-preserving 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

 K=Kc|J0(^ρ)| (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 gyro-averaging operation “breaks” the critical line at the zeros of the zero-order 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

 K≈Kc/|J1(ρi)|^ρ−ρi|, (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 gyro-averaged drift wave map, we introduce the gyroaveraged standard nontwist map (GSNM), which corresponds to a non-monotonic radial electric field. In this case the equilibrium potential is

 ϕ0(x)=α[(xL)−13(xL)3], (24)

where and are dimensional constants. The gyro-average of Eq. (24) gives

 ⟨H0(x)⟩φ=αB0[(xL)(1−ρ22L2)−13(xL)3], (25)

and from Eq. (13) we get the non-monotonic parabolic frequency profile,

 Ω(x)=αB0L⎧⎨⎩⎡⎣1−(ρ√2L)2⎤⎦−(xL)2⎫⎬⎭, (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

 xn+1= xn+2πω0B0AkJ0(kρ)sin(kyn) (27) yn+1= yn+2πω0B0αL⎧⎨⎩⎡⎣1−(ρ√2L)2⎤⎦−(xn+1L)2⎫⎬⎭. (28)

Introducing the dimensionless variables:

 I=−xL,θ=ky2π, (29)

we get the gyro-averaged standard nontwist map (GSNM):

 In+1= In−bJ0(^ρ)sin(2πθn) (30) θn+1= θn+a[(1−¯ρ22)−I2n+1],mod1 (31)

where:

 a= αkω0B0Lb= 2πaAα (32) ¯ρ= ρL^ρ=kρ (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

 bef=bJ0(^ρ). (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

 0= bJ0(^ρ)sin(2πθ∗) (35) m= a[(1−¯ρ22)−I∗2] (36)

where is an integer number. For each , , , and , there are four fixed points:

 P±=(0,±I∗(¯ρ)) (37) Q±=(12,±I∗(¯ρ)) (38)

where

 I∗(¯ρ)=1√2√2(1−ma)−¯ρ2. (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

 R=14[2−Tr(k−1∏i=0J(→xi))] (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

 R(P±)=R(Q∓)=∓πabJ0(^ρ)I∗(¯ρ). (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 del-castillo-1993 (). 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 heteroclinic-type, separatrix reconnection, and the homoclinic-type 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. del-castillo-negrete-etal-1996 (); del-castillo-1993 () and approximate the GSNM in the vicinity of the fixed points by the Hamiltonian

 H(I,θ)=a[I(1−¯ρ22)−I33]−b2πJ0(^ρ)cos(2πθ). (42)

For (configuration I), separatrix reconnection occurs when . If (configuration II), reconnection is observed when . Combining these two conditions, we obtain the reconnection line

 a=34πσ(^ρ,¯ρ)b (43)

with slope:

 σ(^ρ,¯ρ)=|J0(^ρ)|(1−¯ρ22)32 (44)

The reconnection line divides the - parameter space in two regions: one with heteroclinic-type topology () and another with homoclinic-type 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 del-castillo-1993 ()) 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 homoclinic-type 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 heteroclinic-type 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 non-chaotic 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 -degrees-of-freedom Hamiltonian system

 H=H0(I)+ϵV(I,θ,t) (45)

where and are the action-angle variables of the integrable system, , and the parameter controls the strength of the time-periodic perturbation . In this case, the twist (non-degeneracy) condition is

 dΩdI=d2H0dI2≠0, (46)

for all . If there is at least one critical point for which , the non-degeneracy 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 del-castillo-1993 (). 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 (); del-castillo-martinell2013 (). 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 gyro-averaged standard nontwist map

 z10=(14,bJ0(^ρ)2)z20=(34,−bJ0(^ρ)2) (47)
 z11=(a2[1−¯ρ22]+14mod1,0)z21=(a2[1−¯ρ22]+34mod1,0) (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 two-dimensional 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 111Here 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.:

 Rij=Θ(ϵ−∥∥ui−uj∥∥) (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:

 nR=max{n(S(j)R)}Nj=0 (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 zero-order 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 non-chaotic 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. 11a-c. 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 Gyro-averaged quartic nontwist map

In this section, we propose another area-preserving map, the gyro-averaged 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. del-castillo-martinell2013 (). The Hamiltonian system in Ref. del-castillo-martinell2013 () is a drift-wave 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 gyro-averaged drift wave map in Eqs. (11) and (12). The equilibrium potential similar to the one proposed in Ref. del-castillo-martinell2013 ()

 ϕ0(x)=αtanh(xL), (51)

where and are dimensional constants. Applying the gyro-average operation (4) to (51) gives

 ⟨H0(x)⟩φ=α2πB0∫2π0tanh(xL+ρLcosφ)dφ, (52)

and from Eq. (13) we get the nonlinear frequency

 Ω(I)=α2πLB0∫2π0sech2(I+¯ρcosφ)dφ (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 ,

 Ω(I)≈αLB0(1−I2){1−¯ρ2[32(1−I2)−1]}. (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,

 ∂2Ω∂I2∣∣∣I=0,¯ρ=¯ρb=0, (55)

which for Eq. (54) gives

 ¯ρb=1/√2. (56)

Defining , , and , and using Eq. (54), we get from Eqs. (11) and (12), the GQNM

 In+1= In+bJ0(^ρ)sin(2πθn) (57) θn+1= θn+a(1−I2n+1){1−¯ρ2[32(1−I2n+1)−1]},mod1 (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

 0= bJ0(^ρ)sin(2πθ∗) (59) m= a(1−I∗2){1−¯ρ2[32(1−I∗2)−1]},m∈Z. (60)

For , , and , there are four fixed points independently of the value of ,

 P±1=(0,±1),Q±1=(12,±1). (61)

If , there is an additional set of fixed points,

 P±2=(0,±h(¯ρ))Q±2=(12,±h(¯ρ)), (62)

where

 h(¯ρ)=√¯ρ2−23¯ρ2. (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)

 R(P±1) =−R(Q±1)=∓Λ(a,b,^ρ,¯ρ) (64) R(P±2) =−R(Q±2)=±h(¯ρ)Λ(a,b,^ρ,¯ρ) (65)

where

 Λ=πabJ0(^ρ)(1+¯ρ2). (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 non-vanishing 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 -degrees-of-freedom nontwist Hamiltonian systems, the width, , of second-order degenerate resonances scales as

 δΩ∼[ϵ∣∣Ω′′(I0)∣∣]2/3, (67)

where is the amplitude of the perturbation, and is such that

 Ω′(I0)=0Ω′′(I0)≠0, (68)

which according to Eq. (46) corresponds to the point where the twist condition is violated. When the system is slightly perturbed, second-order 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, second-order degenerate resonances might appear near the central shearless curve, and according to Eqs. (67) and (55), the width of these resonances scales as

 lim¯ρ→¯ρbδΩ=0, (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 area-preserving 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 gyro-averaging the electrostatic potential and the maps were constructed from the integration of the resulting gyro-averaged Hamiltonian system.

We constructed and studied three maps: the gyro-averaged standard map (GSM), the gyro-averaged standard nontwist map (GSNM), and the gyro-averaged quartic nontwist map (GQNM). In the GSM, the frequency, i.e. the gyro-averaged zonal flow velocity, has a monotonic radial profile and the map satisfies the twist condition. In the GSNM and GQNM cases the frequency is non-monotonic and as a result the maps are nontwist. Non-monotonic 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 zero-order Bessel function, the breakup diagrams exhibit a highly nontrivial fractal-like 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 period-one orbits. In the GSNM there are four period-one 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 period-one 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/00483-1. 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 UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725. 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. del-Castillo-Negrete, Phys. Plasmas 7, 1702 (2000).
• (5) G. Manfredi and R. O. Dendy, Phys. Rev. Lett., 76, 4360, (1996).
• (6) J.J. Martinell and D. del-Castillo-Negrete, 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. del-Castillo-Negrete and W. Dorland, Phys. Plasmas 16, 102309 (2008).
• (10) D. del-Castillo-Negrete 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. del-Castillo-Negrete 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. del-Castillo-Negrete and P. J. Morrison, Phys. Fluids A 5, 948 (1993)
• (18) F. J. Beron-Vera, 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. del-Castillo-Negrete, 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. de-la-Llave, 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 CLM-PR-12, (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).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters