Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas

A general theoretical framework for investigating nonlinear dynamics of phase space zonal structures is presented in this work. It is then, more specifically, applied to the limit where the nonlinear evolution time scale is smaller or comparable to the wave-particle trapping period. In this limit, both theoretical and numerical simulation studies show that non-adiabatic frequency chirping and phase locking could lead to secular resonant particle transport on meso- or macro-scales. The interplay between mode structures and resonant particles then provides the crucial ingredient to properly understand and analyze the nonlinear dynamics of Alfv\'en wave instabilities excited by non-perturbative energetic particles in burning fusion plasmas. Analogies with autoresonance in nonlinear dynamics and with superradiance in free electron lasers are also briefly discussed.

Alfvénic oscillations can be excited by EPs as well as thermal plasma components [9,10], and, thus, are characterized by a broad spectrum of wavelengths, frequencies and growth rates, which are generally non-trivially interlinked due to the fact that short wavelength kinetic SAW can be excited by resonant mode conversion [11]. Therefore, these drift Alfvén waves (DAW) are expected to play important roles in complex behaviors of burning fusion plasmas [12,13,14,15]; and may have the features of broad band turbulence with O(1) ratio between linear instability growth rate and mode frequency, |γ L /ω 0 | ∼ 1, or nearly periodic (quasi-coherent) nonlinear fluctuations with |γ L /ω 0 | 1 [1]. Meanwhile, such two components of the SAW/DAW spectrum have different effects on EP transports. In this work, the focus is on the nearly periodic (quasi-coherent) nonlinear fluctuations with |γ L /ω 0 | 1. Alfvénic fluctuations in burning plasmas have typically low amplitudes, δ ≡ |δB ⊥ |/B 0 < ∼ 5 × 10 −4 [3], where the ⊥ subscript stands for perpendicular to the equilibrium magnetic field direction b ≡ B 0 /B 0 . Therefore, resonant EPs, more than non resonant ones, are expected to play crucial roles in transport processes [16,17,18,19,20]. Here, in particular, we consider transport processes in two-dimensional (2D) magnetized plasma equilibria with two periodic angle-like coordinates, (θ, ζ), one of which (ζ) defines the equilibrium symmetry. Thus, phase space structures that are of direct interest are those obtained by averaging out dependences on angle-like variables. In fact, in the general procedure used for deriving kinetic equations in weakly-nonideal systems [21,22,23], these structures are those dominating the "principal series" of secular terms obtained by formal perturbation expansion in powers of fluctuating fields; which can be summed up and yields the solution of the corresponding Dyson equation [24,25]. These phase space structures may be generally referred to as "phase space zonal structures" [15,26,27] (PSZS), by analogy with the meso-scale configuration space patterns spontaneously generated by drift-wave turbulence (DWT); i.e., zonal flows and currents/fields or, more generally, zonal structures (ZS) [28,29,30,31]. Here and in the following, meso-scales are intermediate length scales between the SAW/DAW/DWT micro-scales, typically ordered as the perpendicular fluctuation wavelength, λ ⊥ , and the system-size macro-scales, characteristic of equilibrium plasma profiles and macroscopic fluctuations. In magnetized fusion plasmas, at the leading order, ZS are independent of (θ, ζ); due to finite magnetic shear, which suppresses convective cells [32,33]. Hence, ZS do not significantly contribute to cross-field transport [34].
Both ZS and PSZS in fusion plasmas are transverse fluctuations, since they have k ·b = 0 everywhere, with k the wave vector. Meanwhile, they describe the corrugations of equilibrium radial profiles [35,36,37] as consequence of fluctuation induced transport processes, due to emission and reabsorption of (toroidal equilibrium) symmetry breaking perturbations [15]. In general, the time scale of ZS and PSZS evolution is the same as that of the underlying fluctuations, with important consequences on the meso-scale dynamics of SAW/DAW and DWT, which may become nonlocal and be characterized, e.g., by avalanches [38,39,40,41]. For understanding and explaining different behaviors of PSZS, it is important to compare the wave-particle trapping time, τ B , with the characteristic time of nonlinear evolution of PSZS, τ N L . When τ B τ N L , there exists an adiabatic (action) invariant; and phase space density is preserved inside the structure separatrix. When this occurs, phase space holes and clumps [42,43,44,45,46,47,48] are limiting cases of PSZS, and the nonlinear evolution is referred to as "adiabatic". Meanwhile, when τ B ∼ τ N L , no adiabatic (action) invariant exists in the phase space, resonant particle motion can be secular, and the dynamics is "non-adiabatic". Thus, different nonlinear behavior and EP transport are expected, in general, depending on the relative ordering of τ B and τ N L . This also suggests that qualitative and quantitative differences may be expected depending on whether fluctuations are externally controlled, such as in a "driven" system, or are resulting as consequence of plasma instabilities. In fact, while, in the former case, phase space dynamics is typically adiabatic; the evolution, in the instability case, often becomes non-adiabatic such that, at saturation, τ N L > ∼ γ −1 L and τ B > ∼ γ −1 L [27]. In this work, we analyze the case of nonlinear EP distribution functions that dynamically evolve on characteristic "transport" time scales, which are generally of the same order of the nonlinear time scale of the underlying instabilities [27,35,36]; i.e., typically in the non-adiabatic regime. Furthermore, we assume that resonant waveparticle interaction between EPs and Alfvénic fluctuations in burning plasmas is nonperturbative; that is, EPs play significant roles in determining mode frequency and mode structure. Non-perturbative EP response to SAW/DAW instabilities is the crucial element characterizing PSZS nonlinear dynamics discussed in this work. By addressing the general non-adiabatic case, τ N L ∼ γ −1 L ∼ τ B , the present approach can handle both the nonlinear saturation of EP driven Alfvénic instabilities, τ N L ∼ γ −1 L < ∼ τ B , as well as longer time scale nonlinear dynamics, τ N L τ B , such as wave-particle trapping, quasilinear diffusion, and collisional effects. The general theoretical framework is illustrated in Sec. 2, which demonstrates that "radial decoupling", due to resonant EPs exploring finite radial structures of fluctuations in nonuniform plasmas, becomes increasingly more important than and eventually dominates "resonance detuning", due to changing waveparticle phases. Such transition occurs when the size of nonlinear particle orbits is comparable with the radial fluctuation wavelength and corresponds to the onset of nonlocal behaviors; i.e., of meso-scale dynamics. Mode frequency chirping and "phase locking" [16,17,26,27,35,49,50,51,52,53] are shown to be the nonlinear dynamics by which secular resonant particle transport can occur on meso-and macro-scales. In order to illustrate theoretical framework and underlying concepts, Sec. 2 focuses on the obtained results and physics implications, while detailed analytical derivations are given in Appendix A. Furthermore, analogies with "autoresonance" [54,55] in nonlinear systems are pointed out in order to illustrate similarities and differences in the behaviors of "driven" systems and unstable fusion plasmas in the presence of non-perturbative EP free energy sources.
This work aims at presenting fundamental aspects of PSZS nonlinear dynamics to a broad readership; and, thus, presents these issues from different perspectives. For this reason, after discussing the cornerstones of theoretical framework in Sec. 2, we adopt numerical simulations to give examples of how these physics manifest themselves in practice. Thus, meso-scale dynamics of radially localized Energetic Particle Modes (EPMs [56]) are illustrated in Sec. 3, which demonstrates non-perturbative EP behaviors, phase locking and non-adiabatic frequency chirping by numerical simulation results with the extended Hybrid MHD Gyrokinetic Code (XHMGC) [57,58]. Section 4, meanwhile, provides specific derivations and discussions of PSZS nonlinear dynamics within the framework of nonlinear gyrokinetic theory [59,60]. In Sec. 5, the PSZS evolution equations are specialized to the simple case of precessional resonance, for which the nonlinear dynamics is that of 1D non-autonomous nonuniform system. This case allows to derive the general form of the corresponding Dyson equation; and to illustrate the onset of nonlocal behaviors as novel feature in contrast to the case of a 1D uniform plasma. As application of the general theoretical framework of Sec. 5, Sec. 6 analyzes convective amplification of EPM wave packets as avalanches; solitonlike structures in active media, accompanied by secular EP transports on macro-scales. EPM avalanches can be taken as paradigmatic example of "radial decoupling" and nonperturbative EP behaviors. In fact, because of "phase locking", wave-particle resonance is limited only by the finite radial mode structure. At the same time, non-perturbative EP response causes the mode structure to adapt to the modified EP source. Ultimately, the spatiotemporal structure of EPM avalanches is, thus, the self-consistent result of "radial decoupling" and non-perturbative EP response. In Sec. 6, as further evidence of general implications of these physics, analogies with "superradiance" [61] operation regime [62] in free electron lasers (FEL) are also discussed. Section 7 finally gives summarizing remarks and conclusions, while Appendix B and Appendix C are devoted to illustrating detailed expressions and technical derivations, and are proposed to interested readers only.

Theoretical framework
We consider an axisymmetric toroidal magnetized plasma, and use (r, θ, ζ) toroidal flux coordinates with straight magnetic fields lines. Here, r is a radial-like magnetic flux coordinate, ζ is the toroidal symmetry coordinate (cf. Sec. 1), B 0 is represented as in Eq. (A.1), and the magnetic field line pitch q ≡ b · ∇ζ/b · ∇θ = q(r) is a flux function, given in Eq. (A.2). Details about the adopted equilibrium representation and flux coordinates are given in Appendix A. Introducing β = 8πP 0 /B 2 0 as the ratio of kinetic to magnetic pressures, we consider sufficiently low-β plasma equilibria with |∇ ⊥ | |∇ |, such that compressional Alfvén waves are suppressed and magnetic field compression can be solved explicitly via perpendicular pressure balance Thus, δB can be eliminated from governing equations of SAW/DAW and DWT, which are then described in terms of two fluctuating scalar fields; i.e., the scalar potential δφ and the parallel (to b) vector potential δA . As usual, subscript 0 denotes equilibrium quantities and δ is the prefix of fluctuating fields. As anticipated in Sec. 1, we assume that the characteristic nonlinear time, τ N L , satisfies the following ordering [27,63,64] with Ω = eB 0 /(mc) the cyclotron frequency, which generally holds for SAW/DAW excited by EPs in fusion plasmas [15] but also applies to DWT in a variety of cases of practical interest [65,66,67,68]. The ordering of Eq. (2) allows us to address the nonlinear saturation of the considered instabilities, as well as to investigate postsaturation dynamics.

Onset of nonlocal behaviors: resonance detuning and radial decoupling
Based on Eq. (2) and on the low SAW fluctuation amplitudes (cf. Sec. 1), it can be shown that resonant EP motion is slightly modified by fluctuations in one poloidal bounce/transit time. Therefore, fluctuation induced transport is a cumulative effect on resonant EP over many bounce/transit times, which can generally be secular, i.e., not bounded in time. Detailed wave-particle interactions in axisymmetric toroidal magnetized plasmas are analyzed and discussed, under these assumptions, in Appendix A. There, it is shown that the effective action of a generic fluctuation field, f (r, θ, ζ) (time dependences are left implicit), decomposed in Fourier harmonics, can be represented as where ω b is the bounce/transit frequency for magnetically trapped/circulating particles andω d is the toroidal precessional frequency, which are defined in Eqs. (A.5) and (A.11), respectively. Furthermore, τ is a time-like parameter that identifies the particle position along its integrable trajectory in the reference equilibrium, is the bounce harmonic,r is the bounce averaged radial coordinate, defined in Eq. (A.6); and the projection operators P m,n, are defined by Eq. (A.14). Meanwhile, ∆r and Θ m,n, are the nonlinear radial particle displacement and the nonlinear wave-particle phase shift, defined in Eqs. (A. 25) and (A.31), respectively. Mathematically, Eq. (3) is a lifting of f (r, θ, ζ) to the phase space, using the time-like parameter τ to map the effective action of f (r, θ, ζ) as the particle moves along its integrable equilibrium trajectory, identified by its constants of motion. Thus, Eq. (3) reveals the dual nature of fluctuating fields in kinetic descriptions of collisionless plasmas: (i) the field observed in the laboratory frame, which enters in the evolution equation of mode structures; and (ii) the field effectively experienced by the particle in the particle-moving frame, which enters in the description of wave-particle resonances (cf. Appendix A.1). For Θ m,n, = 0 and ∆r = 0, a fluctuation ∝ exp(−iω 0 t) readily yields waveparticle resonance condition in the form of Eqs. (A.16) and (A.17), respectively, for magnetically trapped and circulating particles. Meanwhile, in the presence of fluctuations, Eq. (3) accounts for the cumulative bounce/transit averaged fluctuation effects, discriminating between "resonance detuning", ∝ exp (iΘ m,n, ), and "radial decoupling" ∝ P m,n, • f m,n (r + ∆r). These concepts, introduced theoretically in Ref. [27] and further discussed below and in Appendix A.2, are important for understanding and explaining nonlinear dynamics of SAW/DAW in nonuniform plasmas. In fact, based on the "universal description of a nonlinear resonance" [69] as non-autonomous system with one degree of freedom, which yields a one-dimensional (1D) nonlinear pendulum Hamiltonian also referred to as "standard Hamiltonian" [70], resonance detuning is a general common feature of wave-particle interactions. On the contrary, radial decoupling exists only in nonuniform plasmas; and becomes important when ∆r is comparable with the perpendicular wavelength, λ ⊥ , of perturbation structures ∝ P m,n, • f m,n , which can be estimated as λ ⊥ ∼ ∆ |nq | −1 , ∆ < 1 being a control parameter depending on the considered fluctuation and the nonuniform plasma equilibrium, characterizing the mode transverse scale length [27]. Introducing the "finite interaction length", ∆r L , as the value of ∆r at which wave-particle resonance is lost due to plasma nonuniformity, Eq. (3) shows that a quantitative and qualitative difference in nonlinear dynamics must be expected depending on the relative ordering of ∆r L and λ ⊥ . It is possible to show that ∆r L ∼ 3|γ L /ω 0 |λ −1 n r, with λ n = |nrq | for circulating particles and λ n = 1 for magnetically trapped particles (cf. Appendix A.3, where the corresponding "finite interaction time" is introduced as well). For ∆r L λ ⊥ , the plasma responds like a uniform system, and similarities can be drawn between EP-SAW interactions in burning plasmas of fusion interest and a 1D uniform beam-plasma system near marginality [71]. This limit offers obvious advantages of adopting a simple 1D system for complex dynamics studies [71,72]. However, as demonstrated first in numerical studies of Ref. [73] and further illustrated by recent numerical works [49,74], important differences and novel behaviors of EP-SAW interactions in nonuniform plasmas appear where ∆ < 1 controls the perpendicular fluctuation scale length. Note that Eq. (4) depends on the type of resonance as well as system geometry and nonuniformity through the control parameters λ n and ∆ [26,27]. These numerical simulation studies, illustrate the importance of mode structures in the nonlinear dynamics of SAW/DAW excited by EPs, consistent with Eq. (4) and theoretical predictions [75,76]; and with recent experimental observations in DIII-D, interpreted and explained on the basis of dedicated numerical simulations [77]. Thus, Eq. (4) can be considered as a condition for the onset of nonlocal behaviors in SAW/DAW interactions with EPs; and as indication that the system behaves as a truly 3D nonuniform plasma, despite one isolated resonance can still be described as non-autonomous system with one degree of freedom [15,27]. A systematic investigation of these issues is given in Ref. [53], based on numerical simulation experiments analyzed with phase space diagnostics developed from Hamiltonian mapping techniques.

Phase locking and non-adiabatic phase space dynamics
Equation (4) suggests that onset of nonlocal behaviors and possibly novel features of SAW/DAW-EP interactions with respect to the beam-plasma system could emerge for sufficiently strong EP drive [35]. This poses an interesting issue on whether EP dynamics can be considered perturbatively. For example, the "bump-on-tail paradigm" [27,51,78] (cf. a recent review by Breizman and Sharapov [72]) considers perturbative EP response, which does not modify the thermal plasma dielectric behavior [71]. Thus, in this paradigm, nonlinear dynamics is dominated by wave particle trapping, whose frequency ω B (cf. Appendix A.3) provides the shortest characteristic temporal scale, τ −1 N L . Consequently, there then exists a conserved phase space action connected with wave particle trapping; and nonlinear evolution of EP driven AEs becomes that of phase space holes and clumps, defined as regions with, respectively, lack of density and excess of density w.r.t. the surrounding phase space. These physics have been extensively investigated since the pioneering work by Bernstein, Greene and Kruskal (BGK) [42]; and used for analyzing nonlinear behaviors of 1D uniform Vlasov plasmas [43,44,45,46,47,48], including sources and collisional dissipation [71,79,80]. If the frequency of hole-clump pairs slowly evolves in time [81,82,83,84,85,86], it happens adiabatically at a rate |ω| ω 2 B , set by balancing the rate of energy extraction of the moving holes/clumps in the phase space with the fixed background dissipation ‡. EP transport in velocity space occurs in "buckets" [87], similar to "bucket" EP transport in real space introduced in Ref. [88]. However, unless diffusive transport process are considered due to many overlapping resonances, the EP transport in real space predicted by the "bump-on-tail paradigm" is implicitly local, since it assumes ∆r L λ ⊥ (cf. Sec. 2.1). These transport processes are also similar to "autoresonance" [54,55], by which a nonlinear pendulum can be driven to large amplitude, evolving in time to instantaneously match its frequency with that of an external drive with sufficiently slow downward frequency sweeping.
It was noted by Friedland et al [89] that spontaneous ("autoresonant") nonlinear evolution of phase space structures becomes poorly controllable, when the underlying fluctuations are plasma instabilities rather than being imposed externally. In this respect, non-perturbative EP behaviors become particularly crucial for AEs; when the ‡ Note that the original adiabatic theory of Ref. [84] has been recently extended in Refs. [85,86], but still remains local in the sense discussed in Sec. 2.1 and hereafter. More detailed discussions of this point are given in Refs. [15,27].
following condition is no longer satisfied [63,64]. Here, ∆ω EP is the EP induced complex frequency shift of AEs; and ∆ω SAW is the frequency mismatch between AEs and the closest frequency of the SAW continuous spectrum computed at the peak AE amplitude. When Eq. (5) is progressively broken, AE mode structures are also increasingly affected by EP dynamic response; consistent with theoretical analyses [63,64,75,76] and numerical simulation results [49,73,77,90,91] as well as experimental observations (cf., e.g., Refs. [77,92]). Meanwhile, EP response is always non-perturbative for EPMs [56], excited within the SAW continuous spectrum as discrete fluctuations at the frequency maximizing wave-EP power exchange above the threshold condition set by continuum damping [27,63,64,93]. Clearly, for non-perturbative EP dynamics and significant transports over the domain where fluctuations are localized, the interplay between mode evolution and EP radial redistributions has crucial impact on nonlinear dynamics. Two factors are fundamental for the nonlinear mode evolution: (i) the nonlinear frequency shift determined by EP redistributions, consistent with mode dispersion relation; and (ii) the change in growth rate, taking into account mode structure distortions and EP transports. Both effects are constrained, self-consistently, by the mode dispersive properties. For a given nonlinear frequency shift, there always exists a group of resonant EPs, with given phase w.r.t. the wave, that satisfies |τ N LΘm,n,l |, |τ 2 N LΘ m,n,l | 1 ; i.e., "phase locking"; as shown explicitly for, respectively, circulating and magnetically trapped particles by Eqs. (A.41) and (A.42) (cf. Appendix A.3). These are the resonant EPs that instantaneously maximize wave-particle power exchange [16,17,26,27,35,51,52,53], as they minimize resonance detuning. At the same time, frequency chirping connected with this process is non-adiabatic |ω| < ∼ ω 2 B , as demonstrated by Eqs. (A.37), (A.41) and (A.42). This means that, with non-perturbative EP response, phase locking is signature of non-adiabatic frequency chirping and non-local EP redistributions; and vice versa [27]. Phase locking further extends (but not indefinitely) the finite interaction length and the corresponding finite interaction time (cf. Sec. 2.1 and Appendix A.3) by a factor −1 ω > 1. Here, ω is defined by Eq. (A.45) and depends on the mode dispersive properties. It denotes the effectiveness of phase locking; and, in general, ω < 1 requires the EP dynamics being non-perturbative. Adiabatic chirping (|ω| ω 2 B ) or fixed frequency modes are characterized by ω 1. Taking into account the extended finite interaction length, Eq. (4) can be rewritten as Equation (7) indicates the growth rate threshold for radial mode structures to importantly affect nonlinear dynamics and the transition to non-local EP behaviors. Detailed analysis, would, in general, require a full 3D kinetic description including realistic equilibrium geometry and nonuniformity [26,27]. Assuming typical plasma parameters and the fact that generally ∆ ω < ∼ 10 −1 , Eq. (7) yields the estimates that |γ L /ω 0 | > ∼ 10 −3 for magnetically trapped particles, and |γ L /ω 0 | > ∼ 10 −2 for circulating particles [1,26,27,35,36].
The onset of non-local EP behaviors, above the threshold given by Eq. (7), corresponds to nonlinear meso-scale dynamics involving mode structures [35,73,94,95,96,97]. We discuss these processes in Sec. 3 by numerical simulation experiments, constructed ad hoc to illustrate phase locking and non-adiabatic frequency chirping. However, depending on the considered plasma equilibrium and the corresponding fluctuation spectrum; and for increasing EP drive strength, phase locking and nonlinear meso-scale dynamics can extend to macro-scales. In this way, phase locked EPs could be secularly transported outward by the "mode particle pumping" mechanism, originally introduced for explaining EP loss due to "fishbone" fluctuations [16]. Meanwhile, selfconsistent interplay between nonlinear mode dynamics and EP transports may give rise to EPM "avalanches", by which EPM wave packets are convectively amplified as they propagate radially, until they are quenched by spatial nonuniformity (cf. Sec. 6 and Appendix A.3).

Nonlinear dynamics of localized EPMs
The important roles of plasma non-uniformities and finite radial mode structures for the nonlinear evolution of EP-driven SAW were first discussed by Briguglio et al [73], using hybrid MHD gyrokinetic numerical simulations [98] applied to the analysis of low-n toroidal AE (TAE) and EPM. Further insights into these physics were provided by hybrid MHD gyrokinetic simulations of nonlinear BAE dynamics [49]; employing high-resolution phase space numerical diagnostics techniques based on Hamiltonian mapping [99]. Similar results have also been illustrated by nonlinear gyrokinetic simulations [74].
In order to demonstrate non-perturbative EP behaviors and phase locking more clearly, Briguglio et al [53] have carried out a "numerical simulation experiment" with the XHMGC code [57,58], investigating the nonlinear dynamics of radially localized EPM near the TAE frequency gap. With the aim of controlling the EPM radial localization and inhibiting the onset of EPM convective amplifications (cf. Sec. 6) [35,94,96], Ref. [53] selected a sufficiently flat equilibrium q profile, q(r) = 1.1 + 0.8(r/a) 2 . Adjacent toroidal frequency gaps in the SAW continuous spectrum were, thus, so far apart that nonlocal couplings were minimized, as shown in Fig. 1 (upper panels). For a tokamak equilibrium with major/minor radii ratio R 0 /a = 10, Fig. 1 (upper panels) illustrates the n = 2 SAW continuous spectrum as thin black curve in the background of the intensity contour plots of an n = 2 EPM at different times: ω A t = 1200 (right before mode saturation), ω A t = 1287 (during saturation phase) and ω A t = 1320 (at the end of the EPM burst), with ω A = v A /R 0 computed at the magnetic axis and v A the Alfvén speed. The thermal plasma density, n i (r) = n i0 q 2 (0)/q 2 (r), is chosen such that the SAW frequency gaps are aligned, while the EP guiding-center distribution function (cf. Sec. 4) is a Maxwellian, , ψ(r) is the poloidal magnetic flux function introduced by Eq. (A.1), n EP 0 /n i0 = 1.75 × 10 −3 , v EP /(aΩ EP ) = 0.01 and v EP /v A = 1 (cf. Ref. [53] for further numerical simulation details).
Kinetic Poincaré plots [53] in Fig. 2 show the dynamic evolution of EP PSZS for the same snapshots as in Fig. 1. Particle marker color denotes the initial P φ , normalized to m EP v EP a, being larger (blue) or smaller (red) than the reference resonance value. The wave particle phase, Θ, meanwhile, is normalized to 2π and refers to Θ m,n, (cf. Sec. 2.1) computed for the dominant Fourier and bounce harmonics, (m, n, ) = (3, 2, 1) [53], consistent with the transit resonance condition, Eq. (A.17). The P φ axis is in a.u. and the Θ axis shows the range 0 ≤ Θ 3,2,1 ≤ 4π in order to better visualize EP PSZS. At ω A t = 1200 in Fig. 1(left), the mode is approaching saturation when EPs have completed 1/4 oscillation in the instantaneous fluctuation-induced potential well, as confirmed by the corresponding kinetic Poincaré plot in Fig. 2; and the mode drive is significantly reduced at the original (linear) resonance location. In the further evolution of the mode and after mode saturation, the wave-particle power exchange becomes most negative (damping) at the original resonance location, when resonant EPs have completed 1/2 oscillation in the instantaneous fluctuation-induced potential well, as shown in Fig. 2  for ω A t = 1287. At the same time, due to the distortion of the EP distribution function, steeper gradient regions are formed at the outer limits of the PSZS, as these are the positions where particles tend to accumulate due to radial decoupling and diminishing fluctuation amplitude. This corresponds to an instantaneous strengthening of the mode drive at higher and lower frequencies, with respect to that of the original EPM (damped at this time), which "locally forces" wave-packets within the SAW continuous spectrum in the (r/a, ω/ω A ) plane along the black lines in the background of the top panels of Fig. 1. Thus, the mode structure and frequency tend to split, as shown in both intensity contour plots and in mode structure plots of Fig. 1. Note that both mode saturation as well as mode structure and frequency splitting are non-adiabatic processes, since they occur on a time scale shorter than 2πω −1 B , as demonstrated in Fig. 2. Furthermore, phase locking plays an important role, since the radial dependence of SAW wavepacket frequency and of resonance condition, Eq. (A.17), is essentially the same; and, thus, Eq. (A.41), the phase locking condition, is readily satisfied. Consistent with the analysis of Sec. 2.2, non-adiabatic frequency chirping and phase locking are in one-toone correspondence with non-perturbative EP behaviors. In fact, non-perturbative EP PSZS cause the two distinct dominant Fourier harmonics to shift radially and vary their relative weight; and yield the "rabbit-ear" structure of the (3, 2) mode, visible in the bottom right panel of Fig. 1, that would be otherwise not justifiable.
In the further nonlinear evolution, mode structure and frequency splitting become more evident in Fig. 1 for ω A t = 1320 (right), although the same features are still not visible in Fig. 2. As expected, mode structure splitting appears as double-resonance in the EP PSZS with some delay, as illustrated in Fig. 3 after 3/4 oscillation of resonant EP in the instantaneous fluctuation-induced potential well, due to the time required for EPs to effectively respond to the modified mode structures. After this phase, EP PSZS as well as corresponding mode structures and frequencies tend to merge again, yielding a second EPM burst [53]. As said above, this is due to the ad hoc setup of the present numerical simulation experiment, aimed at inhibiting the onset of EPM convective amplifications (cf. Sec. 6) [35,94,96]. Detailed investigations of these dynamics require an in depth discussion, which will be reported elsewhere. Here, we emphasize again that both intensity contour plots and mode structure distortions are non-perturbative and evidently interlinked self-consistently with the non-adiabatic evolution of EP PSZS (|ω| ∼ ω 2 B ). A crucial role is played by the "radial singular" structures characterizing the SAW continuous spectrum [15,27,63,64]. In fact, while EP response to the regular EPM structure determines the mode frequency [56,63,64] and phase locking (cf. also Sec. 6), the mode structure is peaked ("singular") at the radial position where the SAW continuous spectrum is resonantly excited. Thus, consistent with EP radial redistributions as they are transported nonlinearly, SAW wave-packets readily respond to instantaneous local forcing and maximize wave-particle power exchange. This is consistent with observations from numerical simulation results, reporting that SAW continuous spectrum is crucial in the nonlinear mode dynamics and frequency chirping [73,94,96,97,100,101,102,103,104] (cf. also Appendix A.3).

Governing equations of phase space zonal structures
As prevalent SAW/DAW and DWT fluctuations are in the low-frequency range, |ω 0 | |Ω|, kinetic description is based on the nonlinear gyrokinetic equation. Following Ref.
[59], the perturbed particle distribution function δf can be written as Here, adiabatic and non-adiabatic (δg) particle responses are separated, e −ρ·∇ is the pull-back transformation from guiding-center to particle coordinates [60], ρ ≡ Ω −1 b × v, F 0 is the equilibrium guiding-center particle distribution function, · · · denotes gyrophase averaging, E = v 2 /2 is the energy per unit mass, and µ is the magnetic moment adiabatic invariant (cf. Appendix A). Meanwhile, δg is obtained from the Frieman-Chen nonlinear gyrokinetic equation [59] ∂ ∂t where v d is the magnetic drift velocity κ = b · ∇b is the magnetic field curvature and ∇B 0 κB 0 in the low-β limit. Furthermore, the fluctuation induced particle drift is with δB ⊥g = ∇ × b δA g . Note that Eq. (8) is valid up to O( δ ) and that the relationship between δg and the fluctuating gyrocenter distribution function δF is given by δF = δg + (e/m)∂ EF0 δL g [60]. This approximation is justified within the characteristic (linear/nonlinear) time scale ordering of Eq. (2) (cf. Appendix A.2 for further discussions). Consistent with Sec. 1, the "principal series" of secular terms obtained by formal perturbation expansion in powers of fluctuating fields is dominated by PSZS [21,22,23], which can be summed up and yields the solution of the corresponding Dyson equation [24,25] (cf. Sec. 5). Thus, we formally separate (m, n, ) = (0, 0, 0) terms in Eq. (3) from n = 0 responses. PSZS, denoted hereafter by the subscript z, are then described by the distribution function δf z , decomposed according to Eq. (8); i.e., Here, the double subscript on fields stands for (m, n), J 0 (λ) = e ρ·∇ accounts for finite Larmor radius effects and λ = k ⊥ v ⊥ /Ω. Assuming |k | |k ⊥ |, the PSZS non-adiabatic response can be derived from Eqs. (3) and (10); and is given by [27,35] where ψ is the poloidal magnetic flux, B 0 is represented as in Eq. (A.1), while the first and second terms on the right hand side (RHS) represent, respectively, the effects of ZS (zonal flows/fields) and the corrugation of radial profiles. Note that the formally nonlinear terms ∝ n on the RHS are dominant for n = 0 and |k | |k ⊥ | modes, consistent with the standard gyrokinetic equation orderings [59,60]. However, if the nonlinear interactions between PSZS and ZS themselves are considered, Eq. (14) should be suitably modified and an additional nonlinear term should be added on the RHS (cf. Refs. [51,105]). This is beyond the scope of the present paper, focusing on the nonlinear interaction of PSZS with n = 0 SAW/DAW; and is reported here only for the sake of completeness.
Following the same procedure, the evolution equation for δg n can be derived from Eqs. (3) and (10), yielding On the left hand side (LHS) of Eq. (15), it is interesting to note the effect of ZS in enhancing resonance detuning. Meanwhile, on the RHS, PSZS modify the mode dynamics via corrugation of resonant particle radial profiles. The lack of symmetry between effects of spatial and velocity space gradients of PSZS is only due to having dropped the ∝ δĖ∂ E δg n term on the LHS of Eq. (10), with consistent with the standard gyrokinetic equation orderings [59,60]. Similar to Eq. (14), this term must be restored when considering nonlinear interactions between PSZS and ZS themselves [51,105]. Noting Eq. (A.18), fluctuation induced radial and energy excursion are ordered as ∆r/r ∼ (ω * EP /ω 0 )∆E/E [27,106], with ω * EP the EP diamagnetic frequency and |ω * EP /ω 0 | 1 for EP driven SAW/DAW in fusion plasmas (cf. Appendix A.3). Taking into account the definition of the QF 0 operator in Eq. (16), Eq. (15) can then be rewritten as Here, we have introduced and have consistently neglected O(ω 0 /ω * EP ) terms on the RHS [27]. Using Eq. (19), Eq. (14) can be cast as On the RHS, the dominant contribution in the nonlinear dynamics of EP PSZS is fluctuation induced resonant particle scattering. In Eq. (20), however, the contribution of other collisional [71] and non-collisional [107, 108, 109, 110] scattering processes may be readily included by corresponding scattering operators. Collision induced scattering, introduced in Ref. [71], was the first one to be considered; but microturbulence induced diffusion has been demonstrated to be more relevant in typical ITER conditions [109,110]. Radio-frequency wave induced scattering [107] has also been shown to dominate over collisional scattering in conditions of practical interest in present day experiments [108]. These additional physics affect the nonlinear evolution of F 0EP via self-consistent interplay of SAW/DAW and non-perturbative EP responses on typical time scales longer than those considered here, τ N L ∼ γ −1 L ; and, thus, will be neglected.
Equations (18) and (20) completely describe PSZS nonlinear dynamics, once they are closed by evolution equations for SAW/DAW, δL g n , and ZS, δL g z . They also suggest that nonlinear dynamics of ZS and PSZS may be viewed as generators of neighboring nonlinear equilibria, which generally vary on the same time scale of the underlying fluctuation spectrum [37]. Here, leaving time dependences implicit and using the Poisson summation formula, we follow Ref. [111] and write the spatiotemporal structures of SAW/DAW as [15] δφ n (r, θ, ζ) = 2π with the same decomposition assumed for δA n . Meanwhile, separating parallel mode structure, δφ 0n , and radial envelope, A n (r), we let [15] δφ n (r, ϑ) = A n (r)δφ 0n (r, ϑ) .
Finally, we adopt the theoretical framework of the general fishbone like dispersion relation (GFLDR) [63,64], that allows to write SAW/DAW dynamics equations as which should be intended as the n-th row of a matrix equation for the complex amplitudes (. . . , A n , . . .) involving nonlinear mode couplings [27,36]. Equation (23) is obtained by projection of nonlinear gyrokinetic quasineutrality condition and vorticity equation on the parallel mode structures δφ 0n (r, ϑ); and assumes the eikonal form A n (r) ∼ exp i nq θ k (r)dr [112,113]. Thus, D n (r, θ k , ω) plays the role of a local WKB nonlinear dispersion function, but it can also be written for global fluctuations [63,64]. General expressions of Λ n , δW f n and δW kn are derived in Refs. [63,64] and given in Appendix B without derivation for readers' convenience: Λ n represents a generalized "inertia" accounting for the structures of the SAW continuous spectrum; while δW f n and δW kn can be interpreted as fluid and kinetic "potential energies" [17]. Adopting the time scale ordering of Eq. (2), the spatiotemporal evolution of SAW wave packets can be described expanding the solutions of Eq. (23) about the linear characteristics where D L n is the linearized dispersion function introduced in Eq. (23). Thus, letting becomes [63,111] The S n (r, t) term on the RHS represents a generic source term, including external forcing, S ext n (r, t), and/or nonlinear interactions, denoted by the superscript N L and extracted from the definition of D n in Eq. (23) [27,63,111] (25) also holds as evolution equation for global modes, provided that A n0 is interpreted as mode amplitude at a fixed radial position and D n is suitably redefined as a global dispersion function operator [63] with ω z the ZS complex frequency, Eq. (25) also describes zonal flow dynamics, while the zonal field evolution, assuming massless fluid electron response, is given by [27] with b · ∇δψ ≡ −c −1 ∂ t δA for k = 0 modes. In summary, Eqs. (18) and (20), closed by Eqs. (25) and (27), describe the selfconsistent evolution of ZS and PSZS in the presence of SAW/DAW and/or DWT, based on the time scale ordering of Eq. (2). Despite their relatively simple form, they have been investigated, so far, only in simplified limiting cases; i.e., either neglecting the nonlinear effect of wave-particle resonances [37,65,66,78,114,115,116], or the effects of ZS [35,36,117,118]. In particular, Eqs. (18) and (20) can describe the onset of nonlocal behaviors in SAW/DAW nonlinear dynamics and ensuing EP transport, discussed in Sec. 2.1; and the extension of secular EP transport from meso-to macroscales due to phase locking, described in Sec. 2.2. It was noted in Sec. 2, that the onset of nonlocal behaviors is the signature that, even for an isolated resonance described as a 1D non-autonomous system, the nonuniform plasma dynamics is fully 3D. It is, however, particularly simple and instructive to investigate this problem for precessional resonance only; i.e., Eq. (A.16) with = 0. In such a case, in fact, assuming ω 0 ∼ nω d ω b , EP nonlinear dynamics preserves both µ and J (cf. Appendix A). Thus, the problem is that of a 1D non-autonomous nonuniform system, which can be readily compared with the case of a 1D non-autonomous uniform plasma; e.g., the "bump-on-tail paradigm" (cf. Sec. 2.2).

Dyson equation and the "fishbone paradigm"
In order to simplify the present analysis, we neglect ZS and restrict Eqs. (18) and (20) to precessional resonance with magnetically trapped EPs, while neglecting finite magnetic drift orbit width effects and assuming ω 0 ∼ nω d ω b . In this limit, as anticipated at the end of Sec. 4, we reduce the problem to analyzing a 1D non-autonomous nonuniform system, where onset of nonlocal behaviors is exclusively due to bounce averaged nonlinear EP dynamics. This is the simplest possible case of PSZS evolution, self-consistently interacting with SAW/DAW excited by EPs in fusion plasmas, and has been dubbed as "fishbone paradigm" [27,51,78]. In fact, this is also the case considered originally in the analysis of EP transport and nonlinear burst cycle of resonantly excited internal kink modes [16,17], after their observation as "fishbone" oscillations in the Poloidal Divertor Experiment (PDX) [119]. The peculiar feature of the "fishbone paradigm", emphasized in Refs. [27,51,78], is the self-consistent interplay of mode structures and EP transport, reflecting equilibrium geometry and plasma nonuniformity. In the following, we also show that the "fishbone paradigm" recovers the "bump-on-tail paradigm" in the (local) uniform plasma limit [27,78].

Renormalization of phase space zonal structures
Noting Eq. (3), and due to the frequency ordering ω 0 ∼ nω d ω b as well as to the negligible magnetic drift orbit width limit, we can rewrite Eq. (20) in terms of Thus, recalling that the bounce averaged parallel velocity vanishes for magnetically trapped particles, we have Here, we have used the definition and the ideal MHD condition δE = 0; i.e., δφ = δψ, defined below Eq. (27), to explicitly rewrite the RHS of Eq. (30), formally separating reversible from irreversible processes connected with wave-particle interactions, which dominate the nonlinear dynamics in Eq. (30). In fact, the contribution of reversible processes is O(γ L τ −1 N L /ω 2 0 ) ∼ O(γ 2 L /ω 2 0 ) with respect to irreversible ones, consistent with Eqs. (36) and (45) below. Hence, they will be neglected in the following. Adopting the notation for the Fourier-Laplace transform, Eq. (30) can be solved aŝ Here, we have straightforwardly included the effect of collisions, formally denoted by StF 0 (ω), and of an external source term,Ŝ 0 (ω). FurthermoreF 0 (0) stands for the initial value of F 0 at t = 0, while the repeated n subscript assumes implicit summation. Note, also, that we have explicitly indicated only the dependences on ω (and y , as dummy integration frequency variable) for the sake of notation clarity. Meanwhile, for EP precessional resonance, it is possible to show [27] δK where the subscripts in Q n,yF0 denote mode number and frequency at which the operator defined by Eq. (16) must be evaluated; and, consistent with Eqs. (28) and (29), we have introduced the definition e −inqθ e inqθ ω d δφ n ≡ω dn δφ n .
It is straightforward to verify that Eq. (34) gives back the linear limit forF 0 (ω) = (2πω) −1 iF 0 (0). By substitution of Eq. (34) into Eq. (33), one readily obtainŝ This equation is the analogue of the Dyson equation in quantum field theory (cf., e.g., Ref. [120]) describing EP transport due to emission and reabsorption of (toroidal equilibrium) symmetry breaking perturbations on all possible "loops" [15], schematically shown in Fig. 4. Its solution provides the renormalized expression ofF 0 , and, thereby, of F 0 , taking into account nonuniform toroidal plasmas effects with the addition of sources and collisions. The nonlinear term in Eq. (36) can be further simplified considering that for |ω * EP /ω| 1, typical of SAW/DAW excited by EP in fusion plasmas (cf. Sec. 4).

The uniform plasma limit and the bump-on-tail paradigm
The uniform plasma limit of Eq. (36) illuminates the connection of "fishbone" and "bump-on-tail" paradigms (cf. Secs. 2.2 and 4). It is obtained postulating constant . Schematic view of the Dyson series describing EP transport due to emission and reabsorption of (toroidal equilibrium) symmetry breaking perturbations [15].
δφ n (ω) fluctuations and neglecting spatial dependences ofF 0 by mapping them into corresponding velocity space variations. By direct inspection of Eqs. (15) and (16), and by their comparison with the nonlinear Vlasov equation for a 1D uniform plasma, the configuration (radial) to velocity space mapping is obtained by letting Here, k 0 is the reference wave vector and ω 0 = k 0 v 0 is the resonance condition to be used instead of the precessional resonance of Eq. (A.16) with = 0. Thus, Eq. (38) is completed by the mapping of resonant denominators where L d0 is the characteristic length of variation ofω d § and u ≡ (v − v 0 ). Introducing a source term Q(v) and a simple Krook collision operator, −ν(v)f (v), the Vlasov equation for a 1D uniform plasma equilibrium is written as The configuration to velocity space mapping of Eqs. (38) and (39), thus, allows to rewrite Eqs. (36) and (37) as when expressed for the nonlinear deviation δf 0 (ω) of the particle distribution function from the equilibrium (initial) value F 0 (0) = Q(v)/ν(v). Following Ref. [24,25], it is possible to show that, in the case of many waves characterized by a broad frequency spectrum with overlapping resonances, Eq. (41) reduces to the quasilinear theory of a weakly turbulent plasma [121,122]. In the opposite limit of a narrow frequency spectrum of nearly periodic nonlinear fluctuations (cf. Sec. 5.3), Eq. (41) describes the oscillations of particles that are trapped in the wave, which, however, do not decay in time as expected as consequence of phase mixing. Therefore, the scalar field δφ k 0 oscillations are also predicted to continue indefinitely rather than to fade away, as in actual physical conditions [123], due to the fact that the processes described in Fig. 4 do not account for k 0 -harmonics generation produced by spatial bunching [124]. Equation (41) can also be used to analyze the adiabatic frequency chirping (|ω| ω 2 B ) of phase space hole-clump pairs due to the balance of EP power extraction by PSZS with a fixed background dissipation [81,82,83,84,125]. In particular, two limiting cases have been investigated in the literature: (i) the ω B t 1 limit [81,82,83], where a recursive solution can be adopted near marginal stability; and (ii) the ω B t 1 limit [84,125], where the lowest order solution corresponds to a constant distribution function inside the wave-particle trapping separatrix (hole/clump), which is slowly shifting in velocity space (chirping). The recursive solution of Refs. [81,82,83] which is readily cast as This equation coincides with the evolution equation of PSZS given by Ref. [81], noting that we have introduced ω 4 B ≡ 4(e/m) 2 k 4 0 |δφ k 0 | 2 , in order to preserve the same normalizations of Fourier amplitudes used in the original work. The resulting PSZS dynamics can be shown to admit fixed point solutions, nonlinear oscillations as well as finite time singularities [81,82,83] (cf. Ref. [72] for a recent review). This latter case, obtained for sufficiently low collisionality, corresponds to breaking down of the truncation in the perturbation expansion of the Dyson series [27]; and coincides with the formation of hole/clump pairs, demonstrated in numerical simulations [126,127,128,129]. The analytic description of adiabatic chirping of hole/clump pairs [84,125], valid for ω B t 1 and for isolated PSZS with frequency separation much larger than ω B , shows that the EP distribution function is flattened (in a coarse-grain sense [130]) inside the wave-particle trapping region. Meanwhile, the energy transfer from EP to PSZS during the chirping process is balanced by fixed background dissipation; that is, the condition under which fluctuations are maintained near marginal stability during the adiabatic evolution of hole/clump structures. These features are embedded in the general solution of Eq. (41), which, as noted above, can address wave-particle trapping as the shortest time scale nonlinear dynamics [24,25]. Thus, Eq. (41) can reproduce both ω B t 1 [81,82,83] and ω B t 1 limits [84,125]; and addresses the general ω B t ∼ 1 case, where formation of PSZS is expected [27] (cf. Sec. 3).

Dyson equation for nearly periodic fluctuations
The same results that hold for Eq. (41) with broad and narrow frequency spectrum fluctuations can be obtained for Eq. (36), given the mapping of Eqs. (38) and (39). Thus, one can demonstrate, respectively, quasilinear transport in the radial direction in the case of overlapping resonances, and nonlinear oscillations due to wave-particle trapping in the presence of a nearly periodic large amplitude wave . Further to this, Eq. (36) can address the onset of nonlocal behaviors in EP transport, including effects of equilibrium geometries and plasma nonuniformity discussed in Sec. 2, which are expected to become progressively more important for increasing EP drive and non-perturbative response [27]. The present approach also shows that, even in the perturbative EP local limit and given the mapping of Eqs. (38) and (39), some features remain in the "fishbone" paradigm that make it qualitatively but not quantitatively the same as the "bump-on-tail" paradigm. For example, the energy dependence ofω d andω dn give a different weighting of the particle phase space contributing to the resonant mode drive with respect to the uniform plasma case.
Equation (36), when applied to many modes, provides a framework to explore the transition of EP transports through stochasticity threshold with all the necessary physics ingredients for a realistic comparison with experimental observations and for applications to the dense SAW/DAW fluctuation spectrum typical of fusion plasmas, since it addresses resonance detuning and radial decoupling in wave-particle interactions on the same footing. Equation (36), meanwhile, specialized to nearly periodic fluctuations, can be adopted for nonlinear analyses of "fishbone" oscillations as well as EPMs (cf. Sec. 6).
applies in general to adiabatic as well as non-adiabatic frequency chirping modes. When substituted back into Eq. (36), Eq. (44) yields the following form of the Dyson equation where Q n,ω 0 (τ )F0 can be approximated by Eq. (37). Equation (45) is the renormalized form of PSZS generated by nearly periodic fluctuations, which allows investigating the self-consistent nonlinear evolution of SAW/DAW and EP PSZS when "closed" by Eqs. (18) (dropping the ZS shearing effect ∝ ∂ r δL g z ) and (25). A crucial element, here, is the different argument ofF 0 on the LHS and RHS of Eq. (45), since it bears the information on the self-consistent interplay of nonlinear mode dynamics and EP transport. Assuming |ω| |γ 0 | → 0; i.e., assuming fixed amplitude fluctuations, allows to solve Eq. (45) in the uniform radial limit and to recover nonlinear oscillations due to wave-particle trapping, discussed above in this section. In general, however, we have |ω| ∼ |γ 0 |, and Eq. (45) describes much richer dynamics, among which those to be discussed in Sec. 6.

EPM convective amplification and avalanches
We analyze the nonlinear dynamics of a nearly periodic EPM burst using Eqs. (18), (25) and (45). In particular, we consider a radially localized EPM wave-packet with perpendicular wavelength λ ⊥ L EP , with L EP the characteristic EP pressure scale length. Furthermore, the mode frequency is assumed to be near the lower accumulation point of the TAE frequency gap in the SAW continuum (cf. also Sec. 3). Thus, θ k0 = 0 and ∂D L n /∂θ k0 = 0 in Eq. (25) [63,64,75,76,93]; and Here, the subscript T stands for "toroidal"; and ω /ω u are the lower/upper SAW continuum accumulation point frequency. Using a simple (s, α) model for high aspect ratio (R 0 /a 1) tokamak equilibria with circular shifted magnetic flux surfaces [131], with s ≡ rq /q the magnetic shear and α ≡ −R 0 q 2 β the "ballooning" dimensionless pressure gradient parameter, it is possible to show [63,64,132] δW L f n = |s|π 8 1 + 2κ(s) − α α c ; and 1 2 for |s|, |α| < 1 in Eq. (25), with α c = s 2 /(1+|s|) and κ(s) (1/2)(1+1/|s|) exp(−1/|s|). Meanwhile, consistent with the "fishbone" paradigm (cf. Sec. 5), δW kn can be written as [27,63,64] where λ = µB 0 /E, τ b = 2π/ω b , k ϑ = −nq/r, and we have further assumed that magnetically trapped EPs are characterized by harmonic motion between magnetic mirror points; i.e., they are deeply trapped, such thatω dn nω d in Eq. (35). The linear contribution δW L kn is readily obtained from Eq. (48) by substitution of F 0 (ω) = (2πω) −1 iF 0 (0), as suggested by Eq. (45). In this case, we assume an initial (equilibrium) EP isotropic slowing down distribution function of particles born at energy E F . Thus,F where H denotes the Heaviside step function and the normalization condition is chosen such that the EP energy density is (3/2)P 0EP for E F E c , and EP energy is predominantly transferred to thermal electrons by collisional friction [133] as it occurs for α-particles in fusion plasmas. We then obtain where α EP = −8πR 0 q 2 P 0EP /B 2 0 andω dF = nω d (E = E F /m EP ). Considering a radially localized EP source, α EP = α 0EP exp[−(r − r 0 ) 2 /L 2 EP ] [35]; and using Eqs. (46), (47) and (50), it is possible to show that the linear EPM wave-packet, obtained as solution of Eq. (25), has a characteristic radial envelope width ∼ |L EP /k ϑ | 1/2 , with λ ⊥ ∼ |k ϑ | −1 |L EP /k ϑ | 1/2 L EP , consistent with the initial assumptions. Most important, however, is that resonant EPM-EP interactions play a double role: (i) they give mode drive in excess of the local threshold condition set by SAW continuum damping, ∝ Λ T (ω 0r ) [56]; and (ii) they provide the radial potential well that allows EPM to be radially bounded and excited as absolute instability [64,93]. These properties of resonant EPM-EP interactions are both crucial to explain and understand the nonlinear dynamics of an EPM burst. In fact, they suggest that EPM mode structures may adapt to the nonlinearly modified EP source; and, vice versa, that self-consistent evolution of EP PSZS and EPM mode structures may take the form of avalanches [15,35]. It is, thus, reasonable to seek solutions representing the nearly periodic EPM burst as radially localized, traveling and convectively amplified wave-packets [27].
It is possible to demonstrate that nonlinear dynamics of EP PSZS and EPM are dominated by δW N L kn . In fact, EP effects Λ N L n are negligible due to finite orbit averaging in the kinetic/inertial layer, while δW N L f n accounts for non-resonant EP responses [27,63,64]. Thus, in the present analysis we may assume Λ N L n = δW N L f n = 0. We further set S ext n (r, t) = 0 in Eq. (25); and neglect source and collision terms in Eq. (45), as our aim is to investigate the nonlinear evolution of one single EPM burst excited by the initial unstable distribution given in Eq. (49). By direct substitution into Eq. (48) ofF 0 (ω) from Eq. (45), and keeping the leading order terms in the time scale ordering of Eq. (2), we obtain Note that, here, we have taken into account that nonlinear EPM dynamics occurs on meso-scales, intermediate between λ ⊥ ∼ |k ϑ | −1 and L EP , to have ∂ r commute with quantities that vary on the equilibrium length scale. In Eq. (51), in addition to the resonant denominators, the integrand has two complex-ω dependences that predominantly contribute to the overall integral value: one is connected with the ∝ 1/ω dependence; the other one is due to the nearly singular/peaked structure ofF 0 (ω − 2iγ 0 ) at ω = 2iγ 0 . Due to the ∝ exp(−iωt) dependence of the Laplace integrand, the contribution from ω 0 is exponentially smaller than that from ω 2iγ 0 in the growth and saturation phases of the burst; and is, thus, negligible. This further demonstrates the importance of the self-consistent interplay of nonlinear mode dynamics and EP transport, already noted in Sec. 5.3. Equation (51) also illuminates the competition of phase locking and wave-particle trapping. In fact, noting that the radially localized EPM wave-packet travels at the group velocity v g , there results a relationship between the wave-packet radial position and its time evolving frequency. Thus, focusing on meso-scale radial structures of the integrand in Eq. (51), at the leading order we can rewrite N LΘ m,n, =0 | and ∼ |τ 4 N LΘ m,n, =0 | 2 w.r.t. the first line. These terms are generally important, and, in the "bump-on-tail" paradigm or the uniform plasma limit for SAW/DAW near marginal stability (cf. Sec. 5), account for wave-particle trapping [24,25]. However, when phase locking occurs, Eq. (A.42), second and third line on the RHS of Eq. (52) become negligible, consistent with Eq. (6). Thus, phase locking de facto prevents wave-particle trapping to occur (cf. Appendix A.3). Note also that |τ 2 N LΘ m,n, =0 | ∼ ω 1 by definition, Eq. (A.45). Thus, on time scales shorter than −1 ω τ N L , Eq. (52) with only the first line on the RHS can be used to calculate δW N L kn , which describes convective amplification of a nearly periodic phase locked EPM burst when substituted back into Eq. (25). Phase locking, Eq. (A.42), assumed here as a simplifying Ansatz, can be self-consistently verified a posteriori, once the solution is obtained as shown below. On longer time scales, wave-particle trapping and/or resonance detuning must be accounted for, along with non-resonant EP responses as well as collisions and other relevant phenomena, neglected here but discussed in Ref. [15].
Equation (56) readily recovers the linear limit [64,93]. In the general case [27,63], as suggested by the last term on the RHS, Eq. (56) is a nonlinear Schrödinger equation describing the evolution of the nearly periodic phase locked EPM burst in an "active medium". We can look for solutions in the convectively amplified self-similar form [27] A n0 (r, t) = U (ξ) exp and k n0 denoting the nonlinear radial wave vector. Noting that ∂ rĀn0 = k n0 ∂ ξ U (ξ) exp t γ 0 (t )dt and ∂ tĀn0 = (γ 0 − k n0 v g ∂ ξ )U (ξ) exp t γ 0 (t )dt , there exist two regimes in the nonlinear EPM burst evolution: (i) γ 0 > ∼ |k n0 v g |, where the linear EPM is increasingly modified by the nonlinear interplay with EP transport as the fluctuation intensity increases; and (ii) γ 0 < |k n0 v g |, where the nonlinear EPM-EP interaction dominates the dynamics. In this second case, it is possible to solve Eq. (56) in a simple analytical form. In fact, balancing the nonlinear term with the linear dispersiveness in Eq. (56) yields and withv E×B EP = (−k ϑ c/B 0 ) max[δφ n (r, t)] the radial E ×B velocity at the fluctuation peak, λ g an O(1) control parameter to be determined; and where λ 0 −0.47 + i1.33 corresponds to the ground state of the complex nonlinear oscillator described by Eq. (60), whose numerical solution is shown in Fig. 5, separating amplitude and phase, U (ξ) = W (ξ)e iχ(ξ) . With this solution, Eq. (56) becomes the "nonlinear" EPM dispersion relation where is the linearized "local" EPM dispersion function [56], obtained neglecting the linear dispersiveness ∝ ∂ 2 D L n /∂θ 2 k0 [64]. Equations (59) and (61) describe a one-parameter family, λ g , of possible EPM wavepackets that are convectively amplified as they radially propagate with a group velocity ∝v E×B EP . The dominant nonlinear mode is obtained for i.e., for maximized mode growth and wave-EP power transfer [27]. For typical tokamak plasma parameters, and the simplified model plasma equilibrium adopted here, one generally finds λ g 0.5 ÷ 0.6 + . Thus, the EPM wave-packet propagates radially at a substantial fraction of the peak EP E×B velocity and, as a result, resonant EP motion is predominantly ballistic in the radial direction, similar to the "mode particle pumping" mechanism introduced in Ref. [16] for explaining EP secular loss due to "fishbones". Equation (61) applies at the instantaneous radial position r 1 (t) = r 0 + t 0 v g (t )dt of the EPM wave-packet. The dependence on r 1 , in the present model equilibrium, is explicit in terms ∝ α EP and ∝ω dF ; and implicit in terms depending on ω 0 . Due to the ∝ ln(ω dF /ω 0 − 1) term and ω 0 /ω dF dependences of δW L kn in Eq. (50), the real EPM frequency is given by at the leading oder, so that the phase locking condition, Eqs. (6) and (A.42), is readily satisfied, consistent with the Ansatz and introductory discussions above in this section. Meanwhile, since Imλ 0 > 0, the nonlinear EPM-EP interplay results into a strengthening of the mode drive w.r.t. the linear value ∝ Im δW L kn (ω 0r ) . Thus, the convective amplification of the phase locked EPM is an "avalanche" [35]; and wavepacket amplification can continue until quenched by radial nonuniformities; i.e., the combined effect of decreasing values of linear drive Im δW L kn (ω 0r ) ∝ α EP and increasing values of SAW continuum damping ∝ Λ T (ω 0r ) (cf. Secs. 2.2 and Appendix A.3; and Refs. [15,26,27]). Hence, as consequence of phase locking, meso-scale nonlinear dynamics of the EPM burst as well as EP transport are extended to the macro-scales.
The self-consistent interplay of nonlinear EPM radial structure evolution with EP transport suggests that this is a non-adiabatic "autoresonance" effect (cf. + It is possible to verify that dγ 0 /dλ 2 g > 0 for λ 2 g → 0, and dγ 0 /dλ 2 g < 0 for λ 2 g → ∞. Sec. 2.2), caused by the non-perturbative EP response and the ensuing kinetic EPM dispersiveness. Plasma instability and non-perturbative EP behavior are the crucial differences with respect to adiabatic "autoresonance", where fluctuations are driven/controlled externally and adiabatic frequency sweeping is imposed [54,55]. During the amplification process and non-adiabatic frequency chirping, the relevant nonlinear time scale is τ N L ∼ |∂ −1 t | ∼ |A n0 | −1 , consistent with Eq. (59). Furthermore, since v g = λ gv E×B EP , the EPM intensity is proportional to the square of the traveled distance, similar to the propagation of a FEL short optical pulse in the superradiant regime [62]. Similar to superradiance is also the resonant EPM-EP interaction, which takes place when the wave-EP phase is amplifying due to the phase locking condition. By the time residual resonance detuning changes the wave-EP phase to damping, the resonance condition is lost by radial decoupling. Thus, EPMs can be considered as examples of "autoresonance" and "superradiance" in fusion plasmas [27,134]. Other interesting analogies with neighboring research areas are further discussed in Ref. [27]. For example, Eq. (60) shows interesting analogies with ∂ 2 ξ U = U − 2U 3 yielding U = sech(ξ); i.e., the equation of motion of a nonlinear oscillator in the so-called "Sagdeev potential" V = (−U 2 + U 4 )/2, which can be obtained, e.g., in the analysis of Ion Temperature Gradient turbulence spreading via soliton formation [114], but also from the Gross-Pitaevsky equation, describing the ground state of a quantum system of identical bosons [135,136]; and from the investigation of the envelope of modulated water wave groups [137]. However, Eq. (60) maintains its peculiarities, mostly due to its complex nature coming from the self-consistent interplay of EPM radial structure and EP transport.

Summary and Discussions
Phase space zonal structures play fundamental roles in EP transport induced by SAW/DAW instabilities with |γ L /ω 0 | 1 in burning fusion plasmas. In this work, we have derived general equations that can be used as starting point in analytic as well as numerical investigations of the self-consistent interplay of SAW/DAW instabilities and PSZS in the presence of a non-perturbative EP population in fusion plasmas with realistic equilibrium geometry and plasma nonuniformities. As particular case of the general equations, we have derived a reduced description that may be adopted for nearly periodic fluctuations within the "fishbone" paradigm; and applied it to the analytic study of convective EPM amplification as solitons in active media. Phase locking, in this case, extends the meso-scale dynamics of EPM "avalanches" and causes secular EP redistribution on the macro-scales. These processes have interesting analogies with "autoresonance" in nonlinear dynamics and "superradiance" in FEL operations.
The strength of the instability drive controls the mechanism by which nonlinear saturation is reached in nonuniform plasma equilibria and realistic magnetic field geometries. The reduced 1D uniform plasma description, adopted in the "bump-on-tail" paradigm, applies near marginal stability as long as the nonlinear particle displacement at saturation is significantly smaller than the radial mode width. As a result, mode saturation occurs via wave-particle trapping and resonance detuning. For sufficiently strong EP drive, |γ L /ω 0 | > |γ L /ω 0 | c (|γ L /ω 0 | c < ∼ 10 −2 for typical parameters), the fluctuation induced EP displacement becomes comparable to the radial mode width, and radial decoupling of EPs from mode structures is increasingly more important in the nonlinear saturation dynamics, until it eventually dominates over resonance detuning. This transition is demonstrated in numerical simulation experiments of mesoscale dynamics of EPMs; and is indicative of the plasma behavior as a 3D system, with given equilibrium geometry and plasma nonuniformity. Restricting the analysis to precessional resonance with magnetically trapped particles only, the general case can be reduced to a 1D nonuniform plasma, adopted in the "fishbone" paradigm, which recovers the "bump-on-tail" paradigm in the uniform plasma limit; and allows to understand the onset of non-local behaviors, such as self-consistent interplay of mode structures and EP transport.
The strength of instability drive or, more precisely, the perturbative vs. nonperturbative EP behavior, also controls the frequency dynamics of PSZS. Very near SAW/DAW marginal stability, such that the "bump-on-tail" paradigm can be applied, wave-particle trapping is the fastest nonlinear time scale and, for sufficiently weak collisions, PSZS may evolve into hole/clump pairs with adiabatically sweeping frequency (|ω| ω 2 B ), which are maintained near marginal stability by balancing energy extraction from EP phase space with background dissipation. Meanwhile, for sufficiently strong EP drive that the "fishbone" paradigm must be adopted and EP response is non-perturbative, kinetic SAW/DAW dispersiveness, self-consistently modified by EP transport, causes non-adiabatic frequency evolution (|ω| ∼ ω 2 B ). Phase locked resonant EPs dominate nonlinear dynamics and, by reduction of resonance detuning, nonlinear interplay of secular EP transport with radial mode structures becomes the most important nonlinear physics.
The present theoretical framework, conceived and developed for application to SAW/DAW instabilities in burning fusion plasmas with a non-perturbative EP population, is characterized by the peculiar system geometry and equilibrium nonuniformity. However, it also sheds light on physics aspects that may be of relevance to neighboring fields of research, such as space plasma physics, nonlinear dynamics, condensed matter and accelerator physics. For gyrokinetic theory in axisymmetric toroidal systems, considered in this work, two obvious pairs of action angle coordinates for charged particle motions are (mµ, α), with µ = v 2 ⊥ /(2B 0 )+. . . the magnetic moment (see, e.g., Refs. [60,138]) and α the gyrophase; and (P φ , φ), with P φ the canonical toroidal angular momentum and φ the toroidal angle * . Here, perpendicular (⊥) and parallel ( ) directions are defined with respect to the unit vector b = B 0 /B 0 ; and the axisymmetric equilibrium magnetic field configuration is assumed in the form 2πψ denoting the poloidal magnetic flux within the magnetic surface labeled by ψ. Consistent with the main text, we adopt field aligned toroidal flux coordinates (r, θ, ζ), with r = r(ψ) a radial-like flux variable, the poloidal angle θ is chosen such that the Jacobian J = (∇ψ × ∇θ · ∇ζ) −1 satisfies the condition J B 2 0 being a magnetic flux function [139,140]; and the general toroidal angle ζ ≡ φ − ν(ψ, θ), ν(ψ, θ) being a suitable periodic function of θ [141], such that the safety factor q, is a function of ψ. At the leading order, P φ is given by with Ω = eB 0 /(mc) the cyclotron frequency. The third pair of action angle coordinates is (J, θ c ), with J the "second invariant" and θ c the respective conjugate canonical angle Here, integrals are taken along the particle orbits in the equilibrium magnetic field, dl denoting the arc-length element along it; and we have introduced the unified notation ω b for bounce and transit frequency of trapped and circulating particles, respectively, Given H 0 = m(v 2 + µB 0 ) and the definitions of J and ω b , Therefore, J can be readily computed given the reference magnetic equilibrium J = J(µ, P φ , H 0 ) ↔ H 0 = H 0 (µ, J, P φ ) . * A recent review of coordinates systems and their connection with the description of the guiding-center particle motion is given by Cary and Brizard [138].

Equation (
A.4) yields θ c = ω b τ , with τ a time-like parameter tracking the particle position along the closed poloidal orbit of period 2π/ω b . Thus, for a particle with given constants of motion, wherer andρ c (θ c ) are parameterized by (µ, J, P φ ) and the initial particle radial position ; and "tilde" denotes a generic periodic function in θ c with zero average, which can be computed from the particle equations of motion in the equilibrium B 0 . Similarly, it is readily demonstrated that, for circulating particles (i.e.; particles whose v does not change sign) while, for magnetically trapped particles (i.e.; particles whose v changes sign twice along a closed orbit in the (r, θ) plane), θ =Θ c (θ c ) . Note that, even though q does not depend on θ, it depends on the particle position because of the particle radial displacement. Meanwhile, from Eq. (A.9) and (A.10), one readily derives the following definition for the toroidal precession frequencȳ Appendix A.1. Wave particle resonance condition of a scalar function f (r, θ, ζ), describing a generic fluctuating field for which time dependence has been left implicit, can be projected along the unperturbed particle motion in the equilibrium B 0 using the time-like parameter τ = θ c /ω b ; expressing (r, θ, ζ) The initial particle poloidal angle can be reabsorbed by suitable time shift in the time-like parameter τ ; i.e., in θ c . in terms of the particle position along its unperturbed trajectory. More specifically, we have f (r, θ, ζ) = m,n, ∈Z e i(nω d + ω b )τ P m,n, • f m,n (r) . (A.13) Here, we have introduced the projection operators along the particle motion in the equilibrium B 0 , defined as Furthermore, λ m,n = 1 for magnetically trapped particles, while, for circulating particles, Note that Eq. (A.13) represents the value of f (r, θ, ζ) effectively experienced by the particle along its orbit, parameterized by the actions (µ, J, P φ ), or equivalently by (µ, P φ , H 0 ) (see above), and by the time-like variable τ = θ c /ω b . The integer ∈ Z stands for the "bounce harmonic", similarly to m, n ∈ Z denoting poloidal and toroidal Fourier harmonics. Equation (A.13) is completely defined, together with the functionsΘ c (θ c ), Ξ c (θ c ) andρ c (θ c ) that are obtained from the integration of particle motion, once the reference equilibrium B 0 is assigned. This mode structure decomposition corresponds to a "Lagrangian" description of the fluctuation while moving in the reference frame of the considered particle; i.e., it is a lifting of a generic scalar field to the particle phase space [27]. Meanwhile, Eq. (A.12) is the usual mode structure decomposition of the fluctuating field in the laboratory frame. These two descriptions are the basis of the dual nature of fluctuation structures within a kinetic analysis of wave particle interactions: one with emphasis on the "effective" fluctuation strength experienced by the particles; and another one, reflecting the fields measured by the "observer". Of these two descriptions, the former enters in the resonant particle dynamic response and, thus, is crucially important for wave-particle power exchange, resonant excitation and instability drive; and particle transports [18]. The latter, meanwhile, is what enters in the evolution equation(s) for the mode structures. As noted above, time dependences in f (r, θ, ζ) are left implicit to simplify notation. Once a monochromatic wave f (r, θ, ζ) ∝ exp(−iω 0 t) is assumed, the resonance condition; i.e., the stationarity of wave-particle phase in the particle moving frame where τ ↔ t, is readily derived from Eq. (A.13) and yields for magnetically trapped particles; while, for circulating particles, Here, we remind that ω b stands for both bounce and transit frequencies. We also note thatω d is typically negligible for circulating particles, except for particles that are close to the trapped-to-circulating boundary in the action space.

Appendix A.2. The effect of fluctuations
When considering the effect of fluctuations, for every bounce/transit; i.e., a 2π increment in the angle θ c , a shift will be accumulated in the wave-particle phase ("resonance detuning"), because of radial nonuniformities and the dependence ofω d and ω b on J, P φ that are not conserved anymore in the nonlinear regime (µ is conserved for low frequency fluctuations, analyzed here, for which |ω 0 | |Ω|). Moreover, the wave-particle phase in the mode structure decomposition in action-angle variables will also be shifted, sincẽ Θ c (θ c ),Ξ c (θ c ) andρ c (θ c ) are not simple periodic functions of θ c any longer, as in the given plasma equilibrium. Consistent with the time scale ordering assumed in this work and, in particular, with the condition 1 |ω 0 τ N L | −1 ∼ |γ L /ω 0 | ω ≡ |ω 0 /Ω| (cf. Sec. 2), we assume that for every bounce/transit the effect of the nonlinear dynamics is small compared with the oscillations of the equilibrium particle trajectory. However, the cumulative effect of the nonlinear dynamics on many bounce/transit times can be large and even connected with a secular process; i.e., not bounded in time.
The perturbed equations of motion, which are needed to extend Eq. (A.13) to handle the effect of small but finite amplitude fluctuations, can be readily obtained from the expression of the gyrocenter Hamiltonian up to order ∼ δ ≡ |δB ⊥ /B 0 | 1 [60]; yielding the expressions for δṙ (δψ), δθ and δζ, as well as forḢ 0 ,Ṗ φ andJ. Here, (δr, δθ, δζ) denote the change in particle position due to the small but finite amplitude fluctuations, which are also the cause of δH 0 , δP φ and δJ and finiteḢ 0 , P φ andJ. All these perturbed quantities are considered here to be formally linear in the fluctuation fields. Within this level of approximation of gyrokinetic description of particle motions, the conserved invariants are µ and the extended phase space Hamiltonian K = H 0 + e δφ g − v δA g /c −H (cf., e.g., [70,142]), whose conservation can be readily cast as Furthermore, it is readily shown that, for trapped particles, we have θ =Θ c (θ c ) + ∆θ , (A. 21) while, for circulating particles, while ζ is given by Here, the terms involving ω b and its derivatives with respect to P φ and J are to be considered for circulating particles only. In fact, these terms are originated by the effect of fluctuations on the term ∝qθ c = τ 0 ω bq dτ , due to the unperturbed circulating particle motion, Eq. (A.7), which has no correspondence in the unperturbed trapped particle motion, Eq. (A.8). In particular, terms ∝ ∂ω b /∂P φ , ∂ω b /∂J have the same origin of those included in Eq. (A.22) for circulating particles and excluded in Eq. (A.21) for magnetically trapped particles. Similarly, the term involving ∼ ω b (dq/dr) stems from τ 0 ω b δqdτ for circulating particles only.
With the above results, we can extend the mode structure decomposition in action-angle variables, Eq. (A.13), to include nonlinear particle orbit distortions due to fluctuations:  so that the nonlinear frequency shift ∆ω(t) = ω(t)−ω 0 for a nearly monochromatic wave (cf. Sec. 5.3) is explicitly taken into account, leaving implicit only the time dependence of the reference linear instability. Note that the concept of nearly monochromatic wave still applies for frequency chirping modes as long as |ω| = |∆ω| |γ L ω 0 |; i.e., for both adiabatic (|ω| ω 2 B ) and non-adiabatic (|ω| < ∼ ω 2 B ) frequency sweeping, ω B denoting the wave-particle trapping frequency. Meanwhile, we have, at the lowest order, P m,n, • e ∆r∂r f m,n (r) P m,n, • f m,n (r + ∆r) ; (A. 32) with ∆r being the bounce-averaged nonlinear radial particle displacement, accumulated over many bounces/transits (2π increments in θ c ). For each bounce/transit, the bounceaveraging involved in the definition of P m,n, • f m,n (r), Eq. (A.14), is carried out at the "effective"r location shifted by the accumulated bounce-averaged nonlinear radial particle displacement. Therefore, Eq. (A.28) can be finally cast as Thus, the dominant effect of nonlinear physics on wave-particle resonances enters via the nonlinear phase shift Θ m,n, ("resonance detuning") and/or the P m,n, • f m,n (r + ∆r) functions ("radial decoupling"), which regulate the strength of wave-particle resonances due to the radial mode structure. In both cases, these processes can be understood as cumulative effects of bounce-averaged responses on linear particle motions. This viewpoint is very useful when describing wave-particle interactions by Hamiltonian techniques; e.g., by "kinetic Poincaré plots" introduced recently by White [99]. Meanwhile, a rather detailed analysis, based on implementation of Hamiltonian phase space diagnostics, aimed at investigating the relative importance of "resonance detuning" vs. "radial decoupling" in numerical simulations of nonlinear dynamics of SAW excited by EPs, is given in Ref. [53]. , it is readily recognized that resonance detuning is drastically different for circulating particle and trapped particle resonances. In fact, from EP equations of motion and Eq. (A.18) one generally has ∆r/r ∼ ∆P φ /P φ ∼ (ω * EP /ω 0 )∆H 0 /H 0 [27,106], with ω * EP the EP diamagnetic frequency. Thus, not only radial decoupling but also resonance detuning for EPs is expected to be dominated by the nonlinear radial displacement, since |ω * E /ω 0 | 1 for SAW excited by EPs in fusion plasmas. Let us also consider high-n toroidal mode number modes, as expected in ITER [63,64], and moderate/low bounce harmonics; i.e., = O(1).
Specializing, for simplicity, to the case of a toroidal equilibrium with major radius R 0 and shifted circular magnetic surfaces, we readily obtaiṅ Θ m,n,l n dq dr ω t ∆r − ∆ω , (A. 34) for circulating particles, having denoted the transit frequency as ω t and neglectedω d ; whereas, for magnetically trapped particle resonance, Θ m,n,l n ∂ω d ∂r with λ n = |nrq | for circulating, λ n = 1 for magnetically trapped particles; and δẊ ⊥ the fluctuating E × B guiding-center drift velocity [60,138]. In the derivation of Eq. (A.37), the estimate d tt ∆r ∼ ω 2 B ∆r has been adopted, and resonant particles have been assumed; i.e., Eq. (A.16) for magnetically trapped and Eq. (A.17) for circulating particles, respectively. Given ω B , the separatrix width in particle phase space, ∆r sx , is estimated as ∆r sx |δẊ ⊥ /ω B |. Thus, ∆r sx and ω B are of order ∼ 1/2 δ . Equations (A.34) and (A.35) can be used for estimating the finite interaction length and finite interaction time for resonant particles; i.e., the typical spatial and time scales required for particles to effectively loose the resonance condition in a non-uniform system. For ∆ω = 0; i.e., no frequency chirping, the finite interaction time for resonant particles defines the characteristic nonlinear time scale and is given bẏ Θ m,n, τ N L ∼ 1 ; and τ N L ∼ 1/(3γ L ) . (A.38) Note that the γ −1 L enters via the anticipated saturation process (either resonance detuning or radial decoupling; or both, with varying relative weight); i.e., saturation occurs when τ N L ∼ 1/(3γ L ). Here, the factor 3 is a typical value adopted from analyses of the beam-plasma system [143,144], as well as numerical simulations of Alfvén Eigenmodes (AE) in toroidal plasmas [145]. Meanwhile, neglecting the radial mode structure, the finite interaction length, ∆r L , is obtained from Eq. (A.34) as 3γ L ∼ nq ω∆r L (A. 39) for circulating particles; while, from Eq. (A.35), 3γ L ∼ ω(∆r L /r) (A. 40) for magnetically trapped particles; where characteristic radial scale lengths of ω t ,ω d and ω b are taken to be ∼ r. From these estimates, we readily recognize that the finite interaction length, ∆r L , for circulating resonant particles is typically smaller than the characteristic width, ∼ |nq | −1 [63,64,111], of one poloidal harmonic, δφ m,n ; whereas it may be comparable or larger than ∼ |nq | −1 for magnetically trapped resonant particles if (γ L /ω) > ∼ |3nrq | −1 . For high-n modes as expected in ITER [63,64], this condition is (γ L /ω) > ∼ 10 −2 . Thus, circulating particles transport due to isolated resonances is predicted to be local, causing minor radial profile relaxation, while transport is expected to be mostly diffusive in the presence of many modes [27]. Meanwhile, magnetically trapped EP transports is intrinsically non-local; i.e., characterized by meso-scales, larger than |nq | −1 and shorter than the equilibrium scale length, and convective (ballistic) processes [27]. Different trapped particle vs. circulating particle behaviors are very general, as they are connected with the properties of Eqs. (A.34) and (A. 35), and are expected to play important role whenever resonant particle transport is significant. In fact, both convective and diffusive electron motion have been reported in gyrokinetic numerical simulations of collisionless trapped electron mode turbulence [146]. Ballistic and (super/sub) diffusive EP behaviors have also been demonstrated in the presence of interchange turbulence in simple magnetized toroidal plasmas [147,148]; and, more recently, diffusive circulating particle transport vs. ballistic magnetically trapped particle behaviors induced by electrostatic drift wave turbulence have been discussed in Ref. [149].
Somewhat different behaviors are expected for low-n modes, routinely excited in present day experiments. Here, convective (ballistic) supra-thermal particle transport is observed typically in connection with significant non-adiabatic [|ω| < ∼ ω 2 B ; cf. discussions following Eq. (A.31)] frequency sweeping as, e.g., in NSTX experiments with neutral beam injection [150,151,152]. Secular radial motions of resonant particles were recognized by White et al [16] to be a key mechanism, dubbed "mode particle pumping", for explaining fishbone nonlinear dynamics [17] and related particle losses. This phenomenology is readily explained by Eqs. (A.34) to (A.36), which have the same structure of the original equations derived in Ref. [16]; and predict secular (ballistic) particle motions if theΘ m,n, = 0 resonance condition is preserved as the particles are transported out of the system. In fact, Eqs. (A.39) and (A.40) give "mode particle pumping" by "phase locking" for circulating and magnetically trapped particles when, respectively [16,153] 36) and (A.37), "phase locking" and secular (ballistic) particle transport always imply non-adiabatic frequency sweeping, ∆ω ∼ ω 2 B . Another important feature of the phase locking condition and ∆ṙ given by Eq. (A.36) is that the frequency sweeping rate is proportional to the mode amplitude; i.e., it is fastest when the mode is strongest. This is observed in typical experimental conditions, when resonant particle transport is ballistic (see, e.g., Refs. [3,151]), and in numerical simulations of, e.g., nonlinear EPM evolutions [73,94,95,96,97,100] and, more recently, nonlinear fishbone dynamics [50,154,155].
The effect of phase locking is to extend the finite interaction lengths, given by where numerator and denominator are computed from Eq. (A.31), respectively, at the actual value of the nonlinear frequency shift (frequency chirping) and at ∆ω = 0 (no chirping). The value of ω denotes the effectiveness of phase locking; it depends on the wave dispersive properties [27] and, in general, ω < 1 requires the EP dynamics being non-perturbative. A similar argument applies to the finite interaction time, which is also extended by ∼ −1 ω . Because of the extended interaction length, circulating resonant particles can also be displaced by a significant fraction of ∼ |nq | −1 before resonance detuning sets in. Thus, in the presence of significant non-adiabatic frequency sweeping, circulating EP transports are also expected to occur as convective processes even for long-wavelength modes. This has, indeed, been observed experimentally, e.g., in recent NSTX experiments with neutral beam injection [150,151]; and in numerical simulations [73,94,97,100], where it is shown that the structure of the SAW continuous spectrum is crucial in the nonlinear mode dynamics and frequency chirping [73,94,96,97,100,101,102,103,104].
The secular resonant particle motion, predicted by Eqs. (A.41) and (A.42) is restricted to a group of particles that have similar initial wave-particle phase and can be maintained as long as the mode frequency can be swept preserving phase locking on the one hand and, on the other hand, satisfying the mode dispersion relation. When phase locked resonant particles eventually loose the resonance condition after the finite interaction time/length, they may be substituted by others; and the process can continue until quenched by equilibrium nonuniformity (cf. Secs. 2.2 and 6; and Refs. [15,26,27]). In fact, non-adiabatic frequency sweeping prevents the de facto wave-particle trapping to occur and particle can enter and leave the resonant region during the nonlinear evolution, which maximizes wave-particle power exchange as well as mode growth [27] (cf. Sec. 6).