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

The important role of shear Alfvén waves (SAW) instabilities and their interaction with energetic particles/charged fusion products (EP) in burning fusion plasmas is widely acknowledged (see, e.g., [1–3]), and, historically, it was realized since the pioneering works by Belikov et al [4, 5], Rosenbluth and Rutherford [6] and Mikhailovskii [7, 8].

Alfvénic oscillations can be excited by EPs as well as thermal plasma components [9, 10], and are thus characterized by a broad spectrum of wavelengths, frequencies, and growth rates, which are generally nontrivially 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–15]; and could have the features of broad band turbulence with $\mathcal{O}(1)$ ratio between linear instability growth rate and mode frequency, $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \sim 1$, or nearly periodic quasicoherent nonlinear fluctuations with $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \ll 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 quasicoherent nonlinear fluctuations with $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \ll 1$.

Alfvénic fluctuations in burning plasmas have typically low amplitudes, ${{\epsilon }_{\delta }}\equiv \mid \delta {{{\boldsymbol{B}} }_{\bot }}\mid /{{B}_{0}}\lesssim 5\times {{10}^{-4}}$ [3], where the ⊥ subscript stands for perpendicular to the equilibrium magnetic field direction ${\boldsymbol{b}} \equiv {{{\boldsymbol{B}} }_{0}}/{{B}_{0}}$. Therefore, resonant EPs, more than nonresonant ones, are expected to play crucial roles in transport processes [16–20]. Here, in particular, we consider transport processes in two-dimensional (2D) magnetized plasma equilibria with two periodic angle-like coordinates, $(\theta ,\zeta )$, 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–23], these structures dominate 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 can 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–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, ${{\lambda }_{\bot }}$, 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 $(\theta ,\zeta )$ 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 ${\boldsymbol{k}} \cdot {\boldsymbol{b}} =0$ everywhere, with ${\boldsymbol{k}}$ as the wave vector. Meanwhile, they describe the corrugations of equilibrium radial profiles [35–37] as a 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 could become nonlocal and be characterized, e.g., by avalanches [38–41]. For understanding and explaining different behaviors of PSZS, it is important to compare the wave–particle trapping time, ${{\tau }_{B}}$, with the characteristic time of nonlinear evolution of PSZS, ${{\tau }_{{\rm NL}}}$. When ${{\tau }_{B}}\ll {{\tau }_{{\rm NL}}}$, 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–48] are limiting cases of PSZS, and the nonlinear evolution is referred to as 'adiabatic'. Meanwhile, when ${{\tau }_{B}}\sim {{\tau }_{{\rm NL}}}$, no adiabatic (action) invariant exists in the phase space, resonant particle motion can be secular, and the dynamics is 'nonadiabatic'. Thus, different nonlinear behavior and EP transport are expected, in general, depending on the relative ordering of ${{\tau }_{B}}$ and ${{\tau }_{{\rm NL}}}$. This also suggests that qualitative and quantitative differences might be expected depending on whether fluctuations are externally controlled, such as in a 'driven' system, or are a result of plasma instabilities. In fact, in the former case, phase space dynamics is typically adiabatic, while, in the second case, the evolution is typically nonadiabatic, such that, at saturation, ${{\tau }_{{\rm NL}}}\gtrsim \gamma _{{\rm L}}^{-1}$ and ${{\tau }_{B}}\gtrsim \gamma _{{\rm L}}^{-1}$ [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 nonadiabatic regime. Furthermore, we assume that resonant wave–particle 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. Nonperturbative EP response to SAW/DAW instabilities is the crucial element characterizing PSZS nonlinear dynamics discussed in this work. By addressing the general nonadiabatic case, ${{\tau }_{{\rm NL}}}\sim \gamma _{{\rm L}}^{-1}\sim {{\tau }_{B}}$, the present approach can handle both the nonlinear saturation of EP-driven Alfvénic instabilities, ${{\tau }_{{\rm NL}}}\sim \gamma _{{\rm L}}^{-1}\lesssim {{\tau }_{B}}$, as well as longer time scale nonlinear dynamics, ${{\tau }_{{\rm NL}}}\gg {{\tau }_{B}}$, such as wave–particle trapping, quasi linear diffusion and collisional effects. The general theoretical framework is illustrated in section 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 wave–particle 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–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, section 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 nonperturbative 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 section 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 section 3, which demonstrates nonperturbative EP behaviors, phase locking and nonadiabatic 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 section 5, the PSZS evolution equations are specialized to the simple case of precessional resonance, for which the nonlinear dynamics is that of one-dimensional (1D) nonautonomous nonuniform system. This case allows us to derive the general form of the corresponding Dyson equation; and to illustrate the onset of nonlocal behaviors as a novel feature in contrast to the case of a 1D uniform plasma. As an application of the general theoretical framework of section 5, section 6 analyzes convective amplification of EPM wave packets as avalanches; soliton-like 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, nonperturbative 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 nonperturbative EP response. In section 6, as further evidence of general implications of these physics, analogies with 'super-radiance' [61] operation regime [62] in free-electron lasers (FEL) are also discussed. Section 7 finally gives summarizing remarks and conclusions, while appendices B and C are devoted to illustrating detailed expressions and technical derivations, and are proposed to interested readers only.

We consider an axisymmetric toroidal magnetized plasma, and use $(r,\theta ,\zeta )$ toroidal flux coordinates with straight magnetic field lines. Here, r is a radial-like magnetic flux coordinate, ζ is the toroidal symmetry coordinate (see section 1), ${\boldsymbol{B}}$₀ is represented as in equation (A.1), and the magnetic field line pitch $q\equiv {\boldsymbol{b}} \cdot \nabla \zeta /{\boldsymbol{b}} \cdot \nabla \theta =q(r)$ is a flux function, given in equation (A.2). Details about the adopted equilibrium representation and flux coordinates are given in appendix A. Introducing $\beta =8\pi {{P}_{0}}/B_{0}^{2}$ as the ratio of kinetic to magnetic pressures, we consider sufficiently low-β plasma equilibria with $\mid {{\nabla }_{\bot }}\mid \gg \mid {{\nabla }_{\parallel }}\mid$, such that compressional Alfvén waves are suppressed, and magnetic field compression can be solved explicitly via perpendicular pressure balance

$$\begin{eqnarray}&&{{\nabla }_{\bot }}\left( {{B}_{0}}\delta {{B}_{\parallel }}+4\pi \delta {{P}_{\bot }} \right)\simeq 0.\end{eqnarray} \tag{ 1 }$$

Thus, $\delta {{B}_{\parallel }}$ can be eliminated from the governing equations of SAW/DAW and DWT, which are then described in terms of two fluctuating scalar fields; i.e., the scalar potential $\delta \phi$ and the parallel (to ${\boldsymbol{b}}$) vector potential $\delta {{A}_{\parallel }}$. As usual, subscript 0 denotes equilibrium quantities and δ is the prefix of fluctuating fields.

As anticipated in section 1, we assume that the characteristic nonlinear time, ${{\tau }_{{\rm NL}}}$, satisfies the following ordering [27, 63, 64]

$$\begin{eqnarray}&&1\gg \mid {{\omega }_{0}}{{\tau }_{{\rm NL}}}{{\mid }^{-1}}\sim \mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gg {{\epsilon }_{\omega }}\equiv \mid {{\omega }_{0}}/\Omega \mid ,\end{eqnarray} \tag{ 2 }$$

with $\Omega =e{{B}_{0}}/(mc)$ as 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–68]. The ordering of equation (2) allows us to address the nonlinear saturation of the considered instabilities, as well as to investigate post-saturation dynamics.

2.1. Onset of nonlocal behaviors: resonance detuning and radial decoupling

Based on equation (2) and on the low SAW fluctuation amplitudes (see section 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,\theta ,\zeta )$ (time dependences are left implicit), decomposed in Fourier harmonics, can be represented as

$$\begin{eqnarray}\begin{array}{rcl} f(r,\theta ,\zeta ) & = & \mathop{\sum }\limits_{m,n\in \mathbb{Z}}{{{\rm e}}^{{\rm i}n\zeta -{\rm i}m\theta }}{{f}_{m,n}}(r) \\ {} & \to & \mathop{\sum }\limits_{m,n,\ell \in \mathbb{Z}}{{{\rm e}}^{{\rm i}\left( n{{{\bar{\omega }}}_{d}}+\ell {{\omega }_{b}} \right)\tau +{\rm i}{{\Theta }_{m,n,\ell }}}}{{\mathcal{P}}_{m,n,\ell }}\circ {{f}_{m,n}}(\bar{r}+\Delta r), \\ \end{array}\end{eqnarray} \tag{ 3 }$$

where ${{\omega }_{b}}$ is the bounce/transit frequency for magnetically trapped/circulating particles, and ${{\bar{\omega }}_{d}}$ is the toroidal precessional frequency, which are defined in equations (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, $\bar{r}$ is the bounce-averaged radial coordinate, defined in equation (A.6); and the projection operators ${{\mathcal{P}}_{m,n,\ell }}$ are defined by equation (A.14). Meanwhile, $\Delta r$ and ${{\Theta }_{m,n,\ell }}$ are the nonlinear radial particle displacement and the nonlinear wave-particle phase shift, defined in equations (A.25) and (A.31), respectively. Mathematically, equation (3) is a lifting of $f(r,\theta ,\zeta )$ to the phase space using the time-like parameter τ to map the effective action of $f(r,\theta ,\zeta )$ as the particle moves along its integrable equilibrium trajectory, identified by its constants of motion. Thus, equation (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 (see appendix A.1).

For ${{\Theta }_{m,n,\ell }}=0$ and $\Delta r=0$, a fluctuation $\propto \;{\rm exp} (-{\rm i}{{\omega }_{0}}t)$ readily yields a wave–particle resonance condition in the form of equations (A.16) and (A.17), respectively, for magnetically trapped and circulating particles. Meanwhile, in the presence of fluctuations, equation (3) accounts for the cumulative bounce/transit averaged fluctuation effects, discriminating between 'resonance detuning', $\propto \;{\rm exp} ({\rm i}{{\Theta }_{m,n,\ell }})$, and 'radial decoupling' $\propto \;{{\mathcal{P}}_{m,n,\ell }}\;\circ \;{{f}_{m,n}}(\bar{r}+\Delta r)$. These concepts, introduced theoretically in [27] and further discussed below and in appendix A.2, are important for understanding and explaining the nonlinear dynamics of SAW/DAW in nonuniform plasmas. In fact, based on the 'universal description of a nonlinear resonance' [69] as a nonautonomous system with one degree of freedom, which yields a 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 $\Delta r$ is comparable with the perpendicular wavelength, ${{\lambda }_{\bot }}$, of perturbation structures $\propto \;{{\mathcal{P}}_{m,n,\ell }}\;\circ \;{{f}_{m,n}}$,which can be estimated as ${{\lambda }_{\bot }}\sim {{\epsilon }_{\Delta }}\mid nq^{\prime} {{\mid }^{-1}}$, ${{\epsilon }_{\Delta }}\lt 1$ being a control parameter depending on the considered fluctuation and the nonuniform plasma equilibrium, characterizing the mode perpendicular scale length [27]. Introducing the 'finite interaction length', $\Delta {{r}_{L}}$, as the value of $\Delta r$ at which wave–particle resonance is lost due to plasma nonuniformity, equation (3) shows that a quantitative and qualitative difference in nonlinear dynamics must be expected depending on the relative ordering of $\Delta {{r}_{L}}$ and ${{\lambda }_{\bot }}$. It is possible to show that $\Delta {{r}_{{\rm L}}}\sim 3\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \lambda _{n}^{-1}r$, with ${{\lambda }_{n}}=\mid nrq^{\prime} \mid$ for circulating particles and ${{\lambda }_{n}}=1$ for magnetically trapped particles (see appendix A.3, where the corresponding 'finite interaction time' is introduced as well). For $\Delta {{r}_{L}}\ll {{\lambda }_{\bot }}$, 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 the obvious advantages of adopting a simple 1D system for complex dynamics studies [71, 72]. However, as demonstrated first in numerical studies of [73] and further illustrated by recent numerical works [49, 74], important differences and novel behaviors of EP-SAW interactions in nonuniform plasmas appear for $\Delta {{r}_{L}}\gtrsim {{\lambda }_{\bot }}$; i.e. (see appendix A.3)

$$\begin{eqnarray}&&\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gtrsim \frac{{{\lambda }_{n}}}{3\mid nrq^{\prime} \mid }{{\epsilon }_{\Delta }},\end{eqnarray} \tag{ 4 }$$

where ${{\epsilon }_{\Delta }}\lt 1$ controls the perpendicular fluctuation scale length. Note that equation (4) depends on the type of resonance as well as system geometry and nonuniformity through the control parameters ${{\lambda }_{n}}$ and ${{\epsilon }_{\Delta }}$ [26, 27]. These numerical simulation studies illustrate the importance of mode structures in the nonlinear dynamics of SAW/DAW excited by EPs, consistent with equation (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, equation (4) can be considered as a condition for the onset of nonlocal behaviors in SAW/DAW interactions with EPs; and as an indication that the system behaves as a truly 3D nonuniform plasma, despite one isolated resonance can still be described as nonautonomous system with one degree of freedom [15, 27]. A systematic investigation of these issues is given in [53], based on numerical simulation experiments analyzed with phase space diagnostics developed from Hamiltonian mapping techniques.

2.2. Phase locking and nonadiabatic 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] (see 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 ${{\omega }_{B}}$ (see appendix A.3) provides the shortest characteristic temporal scale, $\tau _{{\rm NL}}^{-1}$. Consequently, there then exists a conserved phase space action connected with wave–particle trapping; and nonlinear evolution of EP-driven Alfvén eigenmodes (AEs) becomes that of phase space holes and clumps, defined as regions with, respectively, a lack of density and excess of density w.r.t. the surrounding phase space. These physics have been extensively investigated since the pioneering work of Bernstein, Greene and Kruskal (BGK) [42]; and used for analyzing nonlinear behaviors of 1D uniform Vlasov plasmas [43–48], including sources and collisional dissipation [71, 79, 80]. If the frequency of hole/clump pairs slowly evolves in time [81–86], it happens adiabatically at a rate $\mid \dot{\omega }\mid \ll \omega _{B}^{2}$, set by balancing the rate of energy extraction of the moving holes/clumps in the phase space with fixed background dissipation⁵ . EP transport in velocity space occurs in 'buckets' [87], similar to 'bucket' EP transport in real space introduced in [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 $\Delta {{r}_{L}}\ll {{\lambda }_{\bot }}$ (see section 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, nonperturbative EP behaviors become particularly crucial for AEs when the following condition

$$\begin{eqnarray}&&\mid \Delta {{\omega }_{{\rm EP}}}\mid \ll \mid \Delta {{\omega }_{{\rm SAW}}}\mid \end{eqnarray} \tag{ 5 }$$

is no longer satisfied [63, 64]. Here, $\Delta {{\omega }_{{\rm EP}}}$ is the EP-induced complex frequency shift of AEs; and $\Delta {{\omega }_{{\rm SAW}}}$ is the frequency mismatch between AEs and the closest frequency of the SAW continuous spectrum computed at the peak AE amplitude. When equation (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 (see, e.g., [77, 92]). Meanwhile, EP response is always nonperturbative 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 nonperturbative 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

$$\begin{eqnarray}&&\mid {{\tau }_{{\rm NL}}}{{\dot{\Theta }}_{m,n,l}}\mid ,\mid \tau _{{\rm NL}}^{2}{{\ddot{\Theta }}_{m,n,l}}\mid \ll 1;\end{eqnarray} \tag{ 6 }$$

i.e., phase locking; as shown explicitly for, respectively, circulating and magnetically trapped particles by equations (A.41) and (A.42) (see appendix A.3). These are the resonant EPs that instantaneously maximize wave–particle power exchange [16, 17, 26, 27, 35, 51–53], as they minimize resonance detuning. At the same time, frequency chirping connected with this process is nonadiabatic $\mid \dot{\omega }\mid \lesssim \omega _{B}^{2}$, as demonstrated by equations (A.37), (A.41) and (A.42). This means that, with nonperturbative EP response, phase locking is signature of nonadiabatic frequency chirping and nonlocal EP redistributions; and vice versa [27]. Phase locking further extends (but not indefinitely) the finite interaction length and the corresponding finite interaction time (see section 2.1 and appendix A.3) by a factor $\epsilon _{{\dot{\omega }}}^{-1}\gt 1$. Here, ${{\epsilon }_{{\dot{\omega }}}}$ is defined by equation (A.45) and depends on the mode-dispersive properties. It denotes the effectiveness of phase locking; and, in general, ${{\epsilon }_{{\dot{\omega }}}}\lt 1$ requires the EP dynamics being nonperturbative. Adiabatic chirping ($\mid \dot{\omega }\mid \ll \omega _{B}^{2}$) or fixed frequency modes are characterized by ${{\epsilon }_{{\dot{\omega }}}}\simeq 1$. Taking into account the extended finite interaction length, equation (4) can be rewritten as

$$\begin{eqnarray}&&\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gtrsim \frac{{{\lambda }_{n}}}{3\mid nrq^{\prime} \mid }{{\epsilon }_{\Delta }}{{\epsilon }_{{\dot{\omega }}}}.\end{eqnarray} \tag{ 7 }$$

Equation (7) indicates the growth rate threshold for radial mode structures to importantly affect nonlinear dynamics and the transition to nonlocal 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 ${{\epsilon }_{\Delta }}{{\epsilon }_{{\dot{\omega }}}}\lesssim {{10}^{-1}}$, equation (7) yields the estimates that $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gtrsim {{10}^{-3}}$ for magnetically trapped particles, and $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gtrsim {{10}^{-2}}$ for circulating particles [1, 26, 27, 35, 36].

The onset of nonlocal EP behaviors, above the threshold given by equation (7), corresponds to nonlinear meso-scale dynamics involving mode structures [35, 73, 94–97]. We discuss these processes in section 3 by numerical simulation experiments, constructed ad hoc to illustrate phase locking and nonadiabatic 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, self-consistent interplay between nonlinear mode dynamics and EP transports could give rise to EPM 'avalanches', by which EPM wave packets are convectively amplified as they propagate radially, until they are quenched by spatial nonuniformity (see section 6 and appendix A.3).

The important roles of plasma nonuniformities 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 beta induced AE (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 nonperturbative 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 (see section 6) [35, 94, 96], reference [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 figure 1 (upper panels). For a tokamak equilibrium with major/minor radii ratio ${{R}_{0}}/a=10$, figure 1 (upper panels) illustrates the n = 2 SAW continuous spectrum as a thin black curve in the background of the intensity contour plots of an n = 2 EPM at different times: ${{\omega }_{{\rm A}}}t=1200$ (right before mode saturation), ${{\omega }_{{\rm A}}}t=1287$ (during saturation phase) and ${{\omega }_{{\rm A}}}t=1320$ (at the end of the EPM burst), with ${{\omega }_{{\rm A}}}={{v}_{{\rm 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 (see section 4) is a Maxwellian, ${{\bar{F}}_{0{\rm EP}}}={{(2\pi {{v}_{{\rm EP}}})}^{-3/2}}{{n}_{E}}(r){\rm exp} [-E/({{m}_{{\rm EP}}}v_{{\rm EP}}^{2})]$, with ${{n}_{{\rm EP}}}(r)={{n}_{{\rm EP}0}}{\rm exp} [-{{\psi }^{2}}(r)/{{\psi }^{2}}(a)]$, $\psi (r)$ is the poloidal magnetic flux function introduced by equation (A.1), ${{n}_{{\rm EP}0}}/{{n}_{i0}}=1.75\times {{10}^{-3}}$, ${{v}_{{\rm EP}}}/(a{{\Omega }_{{\rm EP}}})=0.01$ and ${{v}_{{\rm EP}}}/{{v}_{{\rm A}}}=1$ (see [53] for further numerical simulation details).

Zoom In Zoom Out Reset image size

Kinetic Poincaré plots [53] in figure 2 show the dynamic evolution of EP PSZS for the same snapshots as in figure 1. Particle marker color denotes the initial canonical toroidal angular momentum, ${{P}_{\phi }}$, normalized to ${{m}_{{\rm EP}}}{{v}_{{\rm EP}}}a$, being larger (blue) or smaller (red) than the reference resonance value. The wave–particle phase, Θ, meanwhile, is normalized to $2\pi$ and refers to ${{\Theta }_{m,n,\ell }}$ (see section 2.1) computed for the dominant Fourier and bounce harmonics, $(m,n,\ell )=(3,2,1)$ [53], consistent with the transit resonance condition, equation (A.17). The ${{P}_{\phi }}$ axis is in arbitrary units (a.u.) and the Θ axis shows the range $0\leqslant {{\Theta }_{3,2,1}}\leqslant 4\pi$ in order to better visualize EP PSZS. At ${{\omega }_{{\rm A}}}t=1200$ in figure 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 figure 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 figure 2 for ${{\omega }_{{\rm 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,\omega /{{\omega }_{{\rm A}}})$ plane along the black lines in the background of the top panels of figure 1. Thus, the mode structure and frequency tend to split, as shown in both intensity contour plots and in the mode structure plots of figure 1. Note that both mode saturation, as well as mode structure and frequency splitting, are nonadiabatic processes, since they occur on a time scale shorter than $2\pi \omega _{B}^{-1}$, as demonstrated in figure 2. Furthermore, phase locking plays an important role, since the radial dependence of SAW wave-packet frequency and of the resonance condition [shown in equation (A.17)] are essentially the same; thus, equation (A.41), the phase locking condition, is readily satisfied. Consistent with the analysis of section 2.2, nonadiabatic frequency chirping and phase locking are in one-to-one correspondence with nonperturbative EP behaviors. In fact, nonperturbative 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 figure 1, which would otherwise be unjustifiable.

Zoom In Zoom Out Reset image size

In the further nonlinear evolution, mode structure and frequency splitting become more evident in figure 1 for ${{\omega }_{{\rm A}}}t=1320$ (right), although the same features are still not visible in figure 2. As expected, mode structure splitting appears as double-resonance in the EP PSZS with some delay, as illustrated in figure 3 after $\simeq 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 (see section 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 nonperturbative and evidently interlinked self-consistently with the nonadiabatic evolution of EP PSZS ($\mid \dot{\omega }\mid \sim \omega _{B}^{2}$). 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 (see also section 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 nonlinear mode dynamics and frequency chirping [73, 94, 96, 97, 100–104] (see also appendix A.3).

Zoom In Zoom Out Reset image size

As prevalent SAW/DAW and DWT fluctuations are in the low-frequency range, $\mid {{\omega }_{0}}\mid \ll \mid \Omega \mid$, kinetic description is based on the nonlinear gyrokinetic equation. Following [59], the perturbed particle-distribution function $\delta f$ can be written as

$$\begin{eqnarray}&&\delta f={{e}^{-{\boldsymbol{ \rho }} \cdot \nabla }}\left[ \delta g-\frac{e}{m}\frac{1}{{{B}_{0}}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mu }\left\langle \delta {{L}_{g}} \right\rangle \right]+\frac{e}{m}\left[ \frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}}\delta \phi +\frac{1}{{{B}_{0}}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mu }\delta L \right].\end{eqnarray} \tag{ 8 }$$

Here, adiabatic and nonadiabatic ($\delta g$) particle responses are separated, ${{e}^{-{\boldsymbol{ \rho }} \cdot \nabla }}$ is the pull-back transformation from guiding-center to particle coordinates [60], ${\boldsymbol{ \rho }} \equiv {{\Omega }^{-1}}{\boldsymbol{b}} \times {\boldsymbol{v}}$, $\bar{F}$₀ is the equilibrium guiding-center particle distribution function,

$$\begin{eqnarray}&&\delta {{L}_{g}}=\delta {{\phi }_{g}}-\frac{{{v}_{\parallel }}}{c}\delta {{A}_{\parallel g}}={{e}^{{\boldsymbol{ \rho }} \cdot \nabla }}\delta L={{e}^{{\boldsymbol{ \rho }} }}\cdot \nabla \left( \delta \phi -\frac{{{v}_{\parallel }}}{c}\delta {{A}_{\parallel }} \right),\end{eqnarray} \tag{ 9 }$$

$\left\langle \cdots \right\rangle$ denotes gyrophase averaging, $\mathcal{E}={{v}^{2}}/2$ is the energy per unit mass, and μ is the magnetic moment adiabatic invariant (see appendix A). Meanwhile, $\delta g$ is obtained from the Frieman–Chen nonlinear gyrokinetic equation [59]

$$\begin{eqnarray}\begin{array}{ll} {} & {} & \left( \frac{\partial }{\partial t}+{{v}_{\parallel }}{{\nabla }_{\parallel }}+{{{\boldsymbol{v}} }_{{\rm d}}}\cdot {{\nabla }_{\bot }}+\delta {{\mathop{{\boldsymbol{\bar{X}}} }\limits^{.}}_{\bot }}\cdot {{\nabla }_{\bot }} \right)\delta g \\ {} & {} & \quad =\;-\left( \frac{e}{m}\frac{\partial }{\partial t}\left\langle \delta {{L}_{g}} \right\rangle \frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}}+\frac{c}{{{B}_{0}}}{\boldsymbol{b}} \times \nabla \left\langle \delta {{L}_{g}} \right\rangle \cdot \nabla {{{\bar{F}}}_{0}} \right), \\ \end{array}\end{eqnarray} \tag{ 10 }$$

where ${{{\boldsymbol{v}} }_{{\rm d}}}$ is the magnetic drift velocity

$$\begin{eqnarray}&&{{{\boldsymbol{v}} }_{{\rm d}}}=\frac{{\boldsymbol{b}} }{\Omega }\times \left( \mu \nabla {{B}_{0}}+{\boldsymbol{ \kappa }} v_{\parallel }^{2} \right)\simeq \frac{\left( \mu {{B}_{0}}+v_{\parallel }^{2} \right)}{\Omega }{\boldsymbol{b}} \times {\boldsymbol{ \kappa }} ,\end{eqnarray} \tag{ 11 }$$

${\boldsymbol{ \kappa }} ={\boldsymbol{b}} \cdot \nabla {\boldsymbol{b}}$ is the magnetic field curvature and $\nabla {{B}_{0}}\simeq {\boldsymbol{ \kappa }} {{B}_{0}}$ in the low-β limit. Furthermore, the fluctuation-induced particle drift is

$$\begin{eqnarray}&&\delta {{\mathop{{\boldsymbol{\bar{X}}} }\limits^{.}}_{\bot }}=\frac{c}{{{B}_{0}}}{\boldsymbol{b}} \times \nabla \left\langle \delta {{L}_{g}} \right\rangle +\frac{{{v}_{\parallel }}}{{{B}_{0}}}{\boldsymbol{ \kappa }} \left\langle \delta {{A}_{\parallel g}} \right\rangle =\frac{c}{{{B}_{0}}}{\boldsymbol{b}} \times \nabla \left\langle \delta {{\phi }_{g}} \right\rangle +{{v}_{\parallel }}\frac{\left\langle \delta {{{\boldsymbol{B}} }_{\bot g}} \right\rangle }{{{B}_{0}}},\end{eqnarray} \tag{ 12 }$$

with $\left\langle \delta {{{\boldsymbol{B}} }_{\bot g}} \right\rangle =\nabla \times {\boldsymbol{b}} \left\langle \delta {{A}_{\parallel g}} \right\rangle$. Note that equation (8) is valid up to $\mathcal{O}({{\epsilon }_{\delta }})$ and that the relationship between $\delta g$ and the fluctuating gyrocenter distribution function $\delta \bar{F}$ is given by $\delta \bar{F}=\delta g+(e/m){{\partial }_{\mathcal{E}}}{{\bar{F}}_{0}}\left\langle \delta {{L}_{g}} \right\rangle$ [60]. This approximation is justified within the characteristic (linear/nonlinear) time scale ordering of equation (2) (see appendix A.2 for further discussion).

Consistent with section 1, the 'principal series' of secular terms obtained by formal perturbation expansion in powers of fluctuating fields is dominated by PSZS [21–23], which can be summed up and yields the solution of the corresponding Dyson equation [24, 25] (see section 5). Thus, we formally separate $(m,n,\ell )=(0,0,0)$ terms in equation (3) from $n\ne 0$ responses. PSZS, denoted hereafter by the subscript z, are then described by the distribution function $\delta {{f}_{z}}$, decomposed according to equation (8); i.e.,

$$\begin{eqnarray}\begin{array}{rcl} \delta {{f}_{z}} & = & \mathop{\sum }\limits_{m}\left\{ {{\mathcal{P}}_{m,0,0}}\circ {{\left[ {{J}_{0}}(\lambda )\delta g \right]}_{m,0}} \right\}-{{\left[ {{J}_{0}}(\lambda )\left( \frac{e}{m}\frac{1}{{{B}_{0}}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mu }\left\langle \delta {{L}_{g}} \right\rangle \right) \right]}_{0,0}} \\ {} & {} & +\frac{e}{m}{{\left[ \frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}}\delta \phi +\frac{1}{{{B}_{0}}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mu }\delta L \right]}_{0,0}}. \\ \end{array}\end{eqnarray} \tag{ 13 }$$

Here, the double subscript on fields stands for (m,n), ${{J}_{0}}(\lambda )=\left\langle {{e}^{{\boldsymbol{ \rho }} }}\cdot \nabla \right\rangle$ accounts for finite Larmor radius effects and $\lambda ={{k}_{\bot }}{{v}_{\bot }}/\Omega$. Assuming $\mid {{k}_{\parallel }}\mid \ll \mid {{{\boldsymbol{k}} }_{\bot }}\mid$, the PSZS nonadiabatic response can be derived from equations (3) and (10); and is given by [27, 35]

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial \delta {{g}_{z}}}{\partial t} & = & -{{\mathcal{P}}_{0,0,0}}\circ {{\left( \frac{e}{m}\frac{\partial }{\partial t}{{\left\langle \delta {{L}_{g}} \right\rangle }_{z}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}} \right)}_{0,0}} \\ {} & {} & +{\rm i}\mathop{\sum }\limits_{m}{{\mathcal{P}}_{m,0,0}}\circ \frac{c}{{\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}\mathop{\sum }\limits_{n\ne 0}n{{\left( \delta {{g}_{n}}{{\left\langle \delta {{L}_{g}} \right\rangle }_{-n}} \right)}_{m,0}}, \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where ψ is the poloidal magnetic flux, ${\boldsymbol{B}}$₀ is represented as in equation (A.1), while the first and second terms on the rhs represent, respectively, the effects of ZS (zonal flows/fields) and the corrugation of radial profiles. Note that the formally nonlinear terms $\propto \;n$ on the RHS are dominant for $n\ne 0$ and $\mid {{k}_{\parallel }}\mid \ll \mid {{{\boldsymbol{k}} }_{\bot }}\mid$ modes, consistent with the standard gyrokinetic equation orderings [59, 60]. However, if the nonlinear interactions between PSZS and ZS themselves are considered, equation (14) should be suitably modified and an additional nonlinear term should be added on the rhs (see [51, 105]). This is beyond the scope of the present paper, focusing on the nonlinear interaction of PSZS with $n\ne 0$ SAW/DAW, and is reported here only for the sake of completeness.

Following the same procedure, the evolution equation for $\delta {{g}_{n}}$ can be derived from equations (3) and (10), yielding

$$\begin{eqnarray}&&\begin{array}{l} \left( \frac{\partial }{\partial t}-\frac{{\rm i}nc}{{\rm d}\psi /{\rm d}r}\frac{\partial {{\left\langle \delta {{L}_{g}} \right\rangle }_{z}}}{\partial r}+{{v}_{\parallel }}{{\nabla }_{\parallel }}+{{{\boldsymbol{v}} }_{{\rm d}}}\cdot {{\nabla }_{\bot }} \right)\delta {{g}_{n}} \\ \quad ={\rm i}\frac{e}{m}\left( Q{{{\bar{F}}}_{0}}-\frac{n{{B}_{0}}}{\Omega {\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}{{\mathcal{P}}_{0,0,0}}\delta {{g}_{z}} \right){{\left\langle \circ \delta {{L}_{g}} \right\rangle }_{n}}. \\ \end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&Q{{\bar{F}}_{0}}=\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}}{\rm i}\frac{\partial }{\partial t}+\frac{{\boldsymbol{b}} \times \nabla {{{\bar{F}}}_{0}}}{\Omega }\cdot (-{\rm i}\nabla ).\end{eqnarray} \tag{ 16 }$$

On the lhs of equation (15), note that the effect of ZS enhances resonance detuning. Meanwhile, on the rhs, PSZS modify the mode dynamics via the corrugation of resonant particle radial profiles. The lack of symmetry between the effects of spatial and velocity space gradients of PSZS is only due to having dropped the $\propto \;\delta \dot{\mathcal{E}}{{\partial }_{\mathcal{E}}}\delta {{g}_{n}}$ term on the lhs of equation (10), with

$$\begin{eqnarray}&&\delta \dot{\mathcal{E}}=-\frac{e}{m}{{{\boldsymbol{v}} }_{{\rm d}}}\cdot {{\nabla }_{\bot }}\left\langle \delta {{L}_{g}} \right\rangle ,\end{eqnarray} \tag{ 17 }$$

consistent with the standard gyrokinetic equation orderings [59, 60]. Similar to equation (14), this term must be restored when considering nonlinear interactions between PSZS and ZS [51, 105].

Noting equation (A.18), fluctuation-induced radial and energy excursions are ordered as $\Delta r/r\sim ({{\omega }_{*{\rm EP}}}/{{\omega }_{0}})\Delta \mathcal{E}/\mathcal{E}$ [27, 106], with ${{\omega }_{*{\rm EP}}}$ as the EP diamagnetic frequency and $\mid {{\omega }_{*{\rm EP}}}/{{\omega }_{0}}\mid \gg 1$ for EP-driven SAW/DAW in fusion plasmas (see appendix A.3). Taking into account the definition of the Q $\bar{F}$₀ operator in equation (16), equation (15) can then be rewritten as

$$\begin{eqnarray}&&\left( \frac{\partial }{\partial t}-\frac{{\rm i}nc}{{\rm d}\psi /{\rm d}r}\frac{\partial {{\left\langle \delta {{L}_{g}} \right\rangle }_{z}}}{\partial r}+{{v}_{\parallel }}{{\nabla }_{\parallel }}+{{{\boldsymbol{v}} }_{{\rm d}}}\cdot {{\nabla }_{\bot }} \right)\delta {{g}_{n}}=-\frac{{\rm i}nc}{{\rm d}\psi /{\rm d}r}\frac{\partial {{F}_{0}}}{\partial r}{{\left\langle \delta {{L}_{g}} \right\rangle }_{n}}.\end{eqnarray} \tag{ 18 }$$

Here, we have introduced

$$\begin{eqnarray}&&{{F}_{0}}\equiv {{\bar{F}}_{0}}+{{\mathcal{P}}_{0,0,0}}\circ \delta {{g}_{z}},\end{eqnarray} \tag{ 19 }$$

and have consistently neglected $\mathcal{O}({{\omega }_{0}}/{{\omega }_{*{\rm EP}}})$ terms on the rhs [27]. Using equation (19), equation (14) can be cast as

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial {{F}_{0}}}{\partial t} & = & -\mathcal{P}_{0,0,0}^{2}\circ {{\left( \frac{e}{m}\frac{\partial }{\partial t}{{\left\langle \delta {{L}_{g}} \right\rangle }_{z}}\frac{\partial {{{\bar{F}}}_{0}}}{\partial \mathcal{E}} \right)}_{0,0}} \\ {} & {} & +{\rm i}{{\mathcal{P}}_{0,0,0}}\;\circ \;\mathop{\sum }\limits_{m}{{\mathcal{P}}_{m,0,0}}\circ \frac{c}{{\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}\mathop{\sum }\limits_{n\ne 0}n{{\left( \delta {{g}_{n}}{{\left\langle \delta {{L}_{g}} \right\rangle }_{-n}} \right)}_{m,0}}. \\ \end{array}\end{eqnarray} \tag{ 20 }$$

On the RHS, the dominant contribution in the nonlinear dynamics of EP PSZS is fluctuation-induced resonant particle scattering. In equation (20), however, the contribution of other collisional [71] and noncollisional [107–110] scattering processes could be readily included by corresponding scattering operators. Collision-induced scattering, introduced in [71], was the first type to be considered; however, micro-turbulence-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}_{0{\rm EP}}}$ via self-consistent interplay of SAW/DAW and nonperturbative EP responses on typical time scales longer than those considered here, ${{\tau }_{{\rm NL}}}\sim \gamma _{{\rm L}}^{-1}$; and, thus, they will be neglected.

Equations (18) and (20) completely describe PSZS nonlinear dynamics, once they are closed by evolution equations for SAW/DAW, ${{\left\langle \delta {{L}_{g}} \right\rangle }_{n}}$ and ZS, ${{\left\langle \delta {{L}_{g}} \right\rangle }_{z}}$. They also suggest that the nonlinear dynamics of ZS and PSZS can 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 [111] and write the spatiotemporal structures of SAW/DAW as [15]

$$\begin{eqnarray}\begin{array}{rcl} \delta {{\phi }_{n}}(r,\theta ,\zeta ) & = & 2\pi \mathop{\sum }\limits_{\ell \in \mathbb{Z}}{{{\rm e}}^{{\rm i}n\zeta -{\rm i}nq(\theta -2\pi \ell )}}\delta {{{\hat{\phi }}}_{n}}(r,\theta -2\pi \ell ) \\ {} & = & \mathop{\sum }\limits_{m\in \mathbb{Z}}{{{\rm e}}^{{\rm i}n\zeta -{\rm i}m\theta }}\int {{{\rm e}}^{{\rm i}(m-nq)\vartheta }}\delta {{{\hat{\phi }}}_{n}}(r,\vartheta ){\rm d}\vartheta , \\ \end{array}\end{eqnarray} \tag{ 21 }$$

with the same decomposition assumed for $\delta {{A}_{\parallel n}}$. Meanwhile, separating parallel mode structure, $\delta {{\hat{\phi }}_{0n}}$, and radial envelope, A_n(r), we let [15]

$$\begin{eqnarray}&&\delta {{\hat{\phi }}_{n}}(r,\vartheta )={{A}_{n}}(r)\delta {{\hat{\phi }}_{0n}}(r,\vartheta ).\end{eqnarray} \tag{ 22 }$$

Finally, we adopt the theoretical framework of the general fishbone-like dispersion relation (GFLDR) [63, 64], that allows us to write SAW/DAW dynamics equations as

$$\begin{eqnarray}&&\left[ {\rm i}{{\Lambda }_{n}}-{{\left( \delta {{{\bar{W}}}_{f}}+\delta {{{\bar{W}}}_{k}} \right)}_{n}} \right]{{A}_{n}}(r)={{D}_{n}}(r,{{\theta }_{k}},\omega ){{A}_{n}}(r)=0,\end{eqnarray} \tag{ 23 }$$

which should be intended as the n-th row of a matrix equation for the complex amplitudes $(\ldots ,{{A}_{n}},\ldots )$ involving nonlinear mode couplings [27, 36]. Equation (23) is obtained by the projection of nonlinear gyrokinetic quasineutrality condition and vorticity equation on the parallel mode structures $\delta {{\hat{\phi }}_{0n}}(r,\vartheta )$; and assumes the eikonal form ${{A}_{n}}(r)\sim {\rm exp} \;{\rm i}\int nq^{\prime} {{\theta }_{k}}(r){\rm d}r$ [112, 113]. Thus, ${{D}_{n}}(r,{{\theta }_{k}},\omega )$ plays the role of a local WKB nonlinear dispersion function, but it can also be written for global fluctuations [63, 64]. General expressions of ${{\Lambda }_{n}}$, $\delta {{\bar{W}}_{fn}}$ and $\delta {{\bar{W}}_{kn}}$ are derived in [63, 64] and given in appendix B without derivation for readers' convenience: ${{\Lambda }_{n}}$ represents a generalized 'inertia' accounting for the structures of the SAW continuous spectrum, while $\delta {{\bar{W}}_{fn}}$ and $\delta {{\bar{W}}_{kn}}$ can be interpreted as fluid and kinetic 'potential energies' [17]. Adopting the time-scale ordering of equation (2), the spatiotemporal evolution of SAW wave packets can be described by expanding the solutions of equation (23) about the linear characteristics

$$\begin{eqnarray}&&D_{n}^{{\rm L}}(r,{{\theta }_{k0}}(r),{{\omega }_{0}})=0,\end{eqnarray} \tag{ 24 }$$

where D_n^L is the linearized dispersion function introduced in equation (23). Thus, letting ${{A}_{n}}(r)={\rm exp} (-{\rm i}{{\omega }_{0}}t){{A}_{n0}}(r,t)$, the spatiotemporal evolution equation for ${{A}_{n0}}(r,t)$ becomes [63, 111]

$$\begin{eqnarray}\begin{array}{ll} {} & {} & \frac{\partial D_{n}^{{\rm L}}}{\partial {{\omega }_{0}}}\left( {\rm i}\frac{\partial }{\partial t} \right){{A}_{n0}}+\frac{\partial D_{n}^{{\rm L}}}{\partial {{\theta }_{k0}}}\left( -\frac{{\rm i}}{nq^{\prime} }\frac{\partial }{\partial r}-{{\theta }_{k0}} \right){{A}_{n0}}+\frac{1}{2}\frac{{{\partial }^{2}}D_{n}^{{\rm L}}}{\partial \theta _{k0}^{2}} \\ {} & {} & \quad \times \left[ {{\left( -\frac{{\rm i}}{nq^{\prime} }\frac{\partial }{\partial r}-{{\theta }_{k0}} \right)}^{2}}{{A}_{n0}}-\frac{{\rm i}}{nq^{\prime} }\frac{\partial {{\theta }_{k0}}}{\partial r}{{A}_{n0}} \right]={{S}_{n}}(r,t). \\ \end{array}\end{eqnarray} \tag{ 25 }$$

The ${{S}_{n}}(r,t)$ term on the rhs represents a generic source term, including external forcing, $S_{n}^{{\rm ext}}(r,t)$, and/or nonlinear interactions, denoted by the superscript NL and extracted from the definition of D_n in equation (23) [27, 63, 111]

$$\begin{eqnarray}&&{{S}_{n}}(r,t)=S_{n}^{{\rm ext}}(r,t)-{{{\rm e}}^{{\rm i}{{\omega }_{0}}t}}\left[ {\rm i}\Lambda _{n}^{{\rm NL}}-{{\left( \delta \bar{W}_{f}^{{\rm NL}}+\delta \bar{W}_{k}^{{\rm NL}} \right)}_{n}} \right]{{{\rm e}}^{-{\rm i}{{\omega }_{0}}t}}{{A}_{n0}}.\end{eqnarray} \tag{ 26 }$$

Equation (25) also holds as an evolution equation for global modes, provided that A_n0 is interpreted as the mode amplitude at a fixed radial position and D_n is suitably redefined as a global-dispersion-function operator [63], with ${{\partial }_{r}}{{D}_{n}}={{\partial }_{{{\theta }_{k}}}}{{D}_{n}}=\partial _{{{\theta }_{k}}}^{2}{{D}_{n}}=0$. By letting ${{A}_{n0}}(r,t)\to {{A}_{z}}(r,t)\equiv {{{\rm e}}^{{\rm i}{{\omega }_{z}}t}}\delta {{\phi }_{z}}$, with ${{\omega }_{z}}$ the ZS complex frequency, equation (25) also describes zonal flow dynamics, while the zonal field evolution, assuming a massless fluid electron response, is given by [27]

$$\begin{eqnarray}&&\frac{\partial }{\partial t}\delta {{A}_{\parallel z}}={{\left( \frac{c}{{{B}_{0}}}{\boldsymbol{b}} \times \nabla \delta {{A}_{\parallel }}\cdot \nabla \delta \psi \right)}_{z}},\end{eqnarray} \tag{ 27 }$$

with ${\boldsymbol{b}} \cdot \nabla \delta \psi \equiv -{{c}^{-1}}{{\partial }_{t}}\delta {{A}_{\parallel }}$ for ${{k}_{\parallel }}\ne 0$ modes.

In summary, equations (18) and (20), closed by equations (25) and (27), describe the self-consistent evolution of ZS and PSZS in the presence of SAW/DAW and/or DWT, based on the time-scale ordering of equation (2). Despite their relatively simple form, they have been investigated, so far, only in simplified limiting cases; i.e., neglecting either the nonlinear effect of wave-particle resonances [37, 65, 66, 78, 114–116] or the effects of ZS [35, 36, 117, 118]. In particular, equations (18) and (20) can describe the onset of nonlocal behaviors in SAW/DAW nonlinear dynamics and the ensuing EP transport, discussed in section 2.1, as well as the extension of secular EP transport from meso- to macro-scales due to phase locking, described in section 2.2. It was noted in section 2 that the onset of nonlocal behaviors is a sign that, even for an isolated resonance described as a 1D nonautonomous 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., equation (A.16) with $\ell =0$. In such a case, in fact, assuming ${{\omega }_{0}}\sim n{{\bar{\omega }}_{d}}\ll {{\omega }_{b}}$, EP nonlinear dynamics preserves both μ and J, the second invariant (see appendix A). Thus, the problem is that of a 1D nonautonomous nonuniform system, which can be readily compared with the case of a 1D nonautonomous uniform plasma; e.g., the 'bump-on-tail paradigm' (see section 2.2).

In order to simplify the present analysis, we neglect ZS and restrict equations (18) and (20) to precessional resonance with magnetically trapped EPs, while neglecting finite magnetic drift orbit width effects and assuming ${{\omega }_{0}}\sim n{{\bar{\omega }}_{d}}\ll {{\omega }_{b}}$. In this limit, as anticipated at the end of section 4, we reduce the problem to analyzing a 1D nonautonomous nonuniform system, where the 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 the '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 [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].

5.1. Renormalization of phase space zonal structures

Noting equation (3), and due to the frequency ordering ${{\omega }_{0}}\sim n{{\bar{\omega }}_{d}},\ll {{\omega }_{b}}$ as well as to the negligible magnetic drift orbit width limit, we can rewrite equation (20) in terms of

$$\begin{eqnarray}&&\delta {{\bar{\phi }}_{n}}={{{\rm e}}^{{\rm i}n(\zeta -q\theta )}}\mathop{\sum }\limits_{m}{{\mathcal{P}}_{m,n,0}}\circ \delta {{\phi }_{m,n}}={{{\rm e}}^{-{\rm i}nq\theta }}\overline{{{{\rm e}}^{{\rm i}nq\theta }}\delta {{\phi }_{n}}},\end{eqnarray} \tag{ 28 }$$

with $\overline{(\ldots )}={{(2\pi )}^{-1}}{{\omega }_{b}}\oint (\ldots ){\rm d}\theta /\dot{\theta }$ denoting magnetic drift-bounce averaging (see appendix A), and the corresponding EP response

$$\begin{eqnarray}&&\delta {{\bar{g}}_{n}}={{{\rm e}}^{{\rm i}n(\zeta -q\theta )}}\mathop{\sum }\limits_{m}{{\mathcal{P}}_{m,n,0}}\circ \delta {{g}_{m,n}}={{{\rm e}}^{-{\rm i}nq\theta }}\overline{{{{\rm e}}^{{\rm i}nq\theta }}\delta {{g}_{n}}}.\end{eqnarray} \tag{ 29 }$$

Thus, recalling that the bounce-averaged parallel velocity vanishes for magnetically trapped particles, we have

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial {{F}_{0}}}{\partial t} & = & {\rm i}\mathop{\sum }\limits_{n}\frac{nc}{{\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}\left( \delta {{{\bar{g}}}_{n}}\delta {{{\bar{\phi }}}_{-n}}-\delta {{{\bar{g}}}_{-n}}\delta {{{\bar{\phi }}}_{n}} \right) \\ {} & = & \mathop{\sum }\limits_{n}\frac{nc}{{\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}\left[ \frac{nc}{{\rm d}\psi /{\rm d}r}\frac{\partial {{F}_{0}}}{\partial r}\frac{\partial }{\partial t}{{\mid \frac{\delta {{{\bar{\phi }}}_{n}}}\mid {{{\omega }_{n}}}}^{2}}+{\rm i}\left( \delta {{{\bar{K}}}_{n}}\delta {{{\bar{\phi }}}_{-n}}-\delta {{{\bar{K}}}_{-n}}\delta {{{\bar{\phi }}}_{n}} \right) \right]. \\ \end{array}\end{eqnarray} \tag{ 30 }$$

Here, we have used the definition

$$\begin{eqnarray}&&\delta g\equiv \delta K+{\rm i}(e/m)Q{{\bar{F}}_{0}}\partial _{t}^{-1}\left\langle \delta {{\psi }_{g}} \right\rangle \ \ \end{eqnarray} \tag{ 31 }$$

and the ideal MHD condition $\delta {{E}_{\parallel }}=0$; i.e., $\delta \phi =\delta \psi$, defined below equation (27), to explicitly rewrite the rhs of equation (30), formally separating reversible from irreversible processes connected with wave-particle interactions, which dominate the nonlinear dynamics in equation (30). In fact, the contribution of reversible processes is $\mathcal{O}({{\gamma }_{{\rm L}}}\tau _{{\rm NL}}^{-1}/\omega _{0}^{2})\sim \mathcal{O}(\gamma _{{\rm L}}^{2}/\omega _{0}^{2})$ with respect to irreversible ones, consistent with equations (36) and (45) below. Hence, they will be neglected in the following. Adopting the notation

$$\begin{eqnarray}&&{{\hat{F}}_{0}}(\omega )={{(2\pi )}^{-1}}\int _{0}^{\infty }{{{\rm e}}^{{\rm i}\omega t}}{{F}_{0}}(t){\rm d}t\end{eqnarray} \tag{ 32 }$$

for the Fourier–Laplace transform, equation (30) can be solved as

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{F}}}_{0}}(\omega ) & = & \frac{{\rm i}}{\omega }{\rm St}{{{\hat{F}}}_{0}}(\omega )+\frac{{\rm i}}{\omega }{{{\hat{S}}}_{0}}(\omega )+\frac{{\rm i}}{2\pi \omega }{{{\bar{F}}}_{0}}(0)+\frac{nc}{\omega ({\rm d}\psi /{\rm d}r)} \\ {} & {} & \times \frac{\partial }{\partial r}\int _{-\infty }^{\infty }\left[ \delta {{\widehat{{\bar{\phi }}}}_{n}}(y)\delta {{\widehat{{\bar{K}}}}_{-n}}(\omega -y)-\delta {{\widehat{{\bar{\phi }}}}_{-n}}(y)\delta {{\widehat{{\bar{K}}}}_{n}}(\omega -y) \right]{\rm d}y. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

Here, we have straightforwardly included the effect of collisions, formally denoted by ${\rm St}{{\hat{F}}_{0}}(\omega )$, and of an external source term, ${{\hat{S}}_{0}}(\omega )$. Furthermore, ${{\bar{F}}_{0}}(0)$ stands for the initial value of F₀ 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 a dummy integration frequency variable) for the sake of notation clarity. Meanwhile, for EP precessional resonance, it is possible to show [27]

$$\begin{eqnarray}&&\delta {{\widehat{{\bar{K}}}}_{n}}(\omega )=\frac{e}{m}\int _{-\infty }^{+\infty }\frac{{{{\hat{\omega }}}_{dn}}}{y}\frac{{{Q}_{n,y}}{{{\hat{F}}}_{0}}(\omega -y)}{n{{{\bar{\omega }}}_{d}}-\omega }\delta {{\widehat{{\bar{\phi }}}}_{n}}(y){\rm d}y,\end{eqnarray} \tag{ 34 }$$

where the subscripts in ${{Q}_{n,y}}{{\hat{F}}_{0}}$ denote the mode number and frequency at which the operator defined by equation (16) must be evaluated; and, consistent with equations (28) and (29), we have introduced the definition

$$\begin{eqnarray}&&{{{\rm e}}^{-{\rm i}nq\theta }}\overline{{{{\rm e}}^{{\rm i}nq\theta }}{{\omega }_{d}}\delta {{\phi }_{n}}}\equiv {{\hat{\omega }}_{dn}}\delta {{\widehat{{\bar{\phi }}}}_{n}}.\end{eqnarray} \tag{ 35 }$$

It is straightforward to verify that equation (34) gives back the linear limit for ${{\hat{F}}_{0}}(\omega )={{(2\pi \omega )}^{-1}}{\rm i}{{\bar{F}}_{0}}(0)$. By substitution of equation (34) into equation (33), one readily obtains

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{F}}}_{0}}(\omega ) & = & \frac{{\rm i}}{\omega }{\rm St}{{{\hat{F}}}_{0}}(\omega )+\frac{{\rm i}}{\omega }{{{\hat{S}}}_{0}}(\omega )+\frac{{\rm i}}{2\pi \omega }{{{\bar{F}}}_{0}}(0)+\frac{e}{m}\frac{nc}{\omega ({\rm d}\psi /{\rm d}r)} \\ {} & {} & \times \frac{\partial }{\partial r}\iint _{-\infty }^{\infty }\left[ \delta {{\widehat{{\bar{\phi }}}}_{n}}(y)\frac{{{{\hat{\omega }}}_{d-n}}}{y^{\prime} }\frac{{{Q}_{-n,y^{\prime} }}{{{\hat{F}}}_{0}}(\omega -y-y^{\prime} )}{-n{{{\bar{\omega }}}_{d}}+y-\omega }\delta {{\widehat{{\bar{\phi }}}}_{-n}}(y^{\prime} ) \right. \\ {} & {} & -\left. \delta {{\widehat{{\bar{\phi }}}}_{-n}}(y)\frac{{{{\hat{\omega }}}_{dn}}}{y^{\prime} }\frac{{{Q}_{n,y^{\prime} }}{{{\hat{F}}}_{0}}(\omega -y-y^{\prime} )}{n{{{\bar{\omega }}}_{d}}+y-\omega }\delta {{\widehat{{\bar{\phi }}}}_{n}}(y^{\prime} ) \right]{\rm d}y{\rm d}y^{\prime} . \\ \end{array}\end{eqnarray} \tag{ 36 }$$

This equation is the analogue of the Dyson equation in quantum field theory (see, e.g., [120]) describing EP transport due to emission and reabsorption of (toroidal equilibrium) symmetry breaking perturbations on all possible 'loops' [15], schematically shown in figure 4. Its solution provides the renormalized expression of $\hat{F}$₀, and, thereby, of F₀, taking into account nonuniform toroidal plasma effects with the addition of sources and collisions. The nonlinear term in equation (36) can be further simplified when considering that

$$\begin{eqnarray}&&{{Q}_{n,\omega }}{{\hat{F}}_{0}}\simeq -\frac{n{{B}_{0}}}{\Omega {\rm d}\psi /{\rm d}r}\frac{\partial {{{\hat{F}}}_{0}}}{\partial r}\end{eqnarray} \tag{ 37 }$$

for $\mid {{\omega }_{*{\rm EP}}}/\omega \mid \gg 1$, typical of SAW/DAW excited by EP in fusion plasmas (see section 4).

Zoom In Zoom Out Reset image size

5.2. The uniform plasma limit and the bump-on-tail paradigm

The uniform plasma limit of equation (36) illuminates the connection of the 'fishbone' and 'bump-on-tail' paradigms (see sections 2.2 and 4). It is obtained by postulating constant $\delta {{\widehat{{\bar{\phi }}}}_{n}}(\omega )$ fluctuations and neglecting spatial dependences of $\hat{F}$₀ by mapping them into corresponding velocity space variations. By direct inspection of equations (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

$$\begin{eqnarray}&&-\frac{n{{B}_{0}}}{{\rm d}\psi /{\rm d}r}\frac{\partial }{\partial r}\leftrightarrow {{k}_{0}}\frac{\partial }{\partial v}.\end{eqnarray} \tag{ 38 }$$

Here, k₀ is the reference wave vector and ${{\omega }_{0}}={{k}_{0}}{{v}_{0}}$ is the resonance condition to be used instead of the precessional resonance of equation (A.16) with $\ell =0$. Thus, equation (38) is completed by the mapping of resonant denominators

$$\begin{eqnarray}&&n{{\bar{\omega }}_{d}}-{{\omega }_{0}}\simeq n{{\bar{\omega }}_{d0}}\frac{(r-{{r}_{0}})}{{{L}_{d0}}}\leftrightarrow {{k}_{0}}(v-{{v}_{0}})={{k}_{0}}u,\end{eqnarray} \tag{ 39 }$$

where L_d0 is the characteristic length of variation of ${{\bar{\omega }}_{d}}$ ⁶ and $u\equiv (v-{{v}_{0}})$.

Introducing a source term Q(v) and a simple Krook collision operator, $-\nu (v)f(v)$, the Vlasov equation for a 1D uniform plasma equilibrium is written as

$$\begin{eqnarray}&&{{d}_{t}}f(v)=Q(v)-\nu (v)f(v).\end{eqnarray} \tag{ 40 }$$

The configuration to velocity space mapping of equations (38) and (39) thus allows us to rewrite equations (36) and (37) as

$$\begin{eqnarray}\begin{array}{rcl} (-{\rm i}\omega +\nu )\delta {{{\hat{f}}}_{0}}(\omega ) & = & {\rm i}\frac{{{e}^{2}}k_{0}^{2}}{{{m}^{2}}}\frac{\partial }{\partial u}\iint _{-\infty }^{\infty }\left[ \delta {{{\hat{\phi }}}_{{{k}_{0}}}}(y)\frac{-{{\partial }_{u}}{{{\hat{F}}}_{0}}(\omega -y-y^{\prime} )}{y-{{k}_{0}}u-\omega -{\rm i}\nu }\delta {{{\hat{\phi }}}_{-{{k}_{0}}}}(y^{\prime} ) \right. \\ {} & {} & -\left. \delta {{{\hat{\phi }}}_{-{{k}_{0}}}}(y)\frac{{{\partial }_{u}}{{{\hat{F}}}_{0}}(\omega -y-y^{\prime} )}{y+{{k}_{0}}u-\omega -{\rm i}\nu }\delta {{{\hat{\phi }}}_{{{k}_{0}}}}(y^{\prime} ) \right]{\rm d}y{\rm d}y^{\prime} , \\ \end{array}\end{eqnarray} \tag{ 41 }$$

when expressed for the nonlinear deviation $\delta {{\hat{f}}_{0}}(\omega )$ of the particle-distribution function from the equilibrium (initial) value ${{F}_{0}}(0)=Q(v)/\nu (v)$. Following [24, 25], it is possible to show that, in the case of many waves characterized by a broad frequency spectrum with overlapping resonances, equation (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 (see section 5.3), equation (41) describes the oscillations of particles that are trapped in the wave, which, however, do not decay in time as expected as a consequence of phase mixing. Therefore, in the scalar field $\delta {{\phi }_{{{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 figure 4 do not account for k₀-harmonics generation produced by spatial bunching [124].

Equation (41) can also be used to analyze the adiabatic frequency chirping ($\mid \dot{\omega }\mid \ll \omega _{B}^{2}$) of phase space hole/clump pairs due to the balance of EP power extraction by PSZS with a fixed background dissipation [81–84, 125]. In particular, two limiting cases have been investigated in the literature: (i) the ${{\omega }_{B}}t\ll 1$ limit [81–83], where a recursive solution can be adopted near marginal stability; and (ii) the ${{\omega }_{B}}t\gg 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 [81–83], obtained for ${{\omega }_{B}}t\ll 1$, corresponds to taking ${{\hat{F}}_{0}}(\omega -y-y^{\prime} )={\rm i}{{(2\pi )}^{-1}}{{F}_{0}}(0){{(\omega -y-y^{\prime} )}^{-1}}$ in equation (41); i.e., to considering only the first loop in the Dyson series, schematically shown in figure 4. Moving to the t-representation, the recursive solution of equation (41) is then obtained as

$$\begin{eqnarray}\begin{array}{rcl} \left( \frac{\partial }{\partial t}+\nu \right)\delta {{f}_{0}} & = & {\rm i}\frac{{{e}^{2}}k_{0}^{2}}{{{m}^{2}}}\frac{\partial }{\partial u}\iint _{-\infty }^{\infty }{{{\rm e}}^{-{\rm i}(y+{{y}^{\prime }})t}}\left[ \delta {{{\hat{\phi }}}_{{{k}_{0}}}}(y)\frac{{{\partial }_{u}}{{F}_{0}}(0)}{y^{\prime} +{{k}_{0}}u+{\rm i}\nu }\delta {{{\hat{\phi }}}_{-{{k}_{0}}}}(y^{\prime} ) \right. \\ {} & {} & +\left. \delta {{{\hat{\phi }}}_{-{{k}_{0}}}}(y)\frac{{{\partial }_{u}}{{F}_{0}}(0)}{y^{\prime} -{{k}_{0}}u+i\nu }\delta {{{\hat{\phi }}}_{{{k}_{0}}}}(y^{\prime} ) \right]{\rm d}y{\rm d}y^{\prime} , \\ \end{array}\end{eqnarray} \tag{ 42 }$$

which is readily cast as

$$\begin{eqnarray}&&\left( \frac{\partial }{\partial t}+\nu \right)\delta {{f}_{0}}=\frac{\omega _{B}^{2}(t)}{4k_{0}^{2}}\frac{\partial }{\partial u}\int _{0}^{t}\left[ {{{\rm e}}^{-(\nu +{\rm i}{{k}_{0}}u)(t-{{t}^{\prime }})}}+c.c. \right]\omega _{B}^{2}(t^{\prime} )\frac{\partial {{F}_{0}}(0)}{\partial u}{\rm d}t^{\prime} .\end{eqnarray} \tag{ 43 }$$

This equation coincides with the evolution equation of PSZS given by [81], noting that we have introduced $\omega _{B}^{4}\equiv 4{{(e/m)}^{2}}k_{0}^{4}\mid \delta {{\phi }_{{{k}_{0}}}}{{\mid }^{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–83] (see [72] for a recent review). This latter case, obtained for sufficiently low collisionality, corresponds to the 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–129]. The analytic description of the adiabatic chirping of hole/clump pairs [84, 125], valid for ${{\omega }_{B}}t\gg 1$ and for isolated PSZS with frequency separation much larger than ${{\omega }_{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 equation (41), which, as noted above, can address wave-particle trapping as the shortest time-scale nonlinear dynamics [24, 25]. Thus, equation (41) can reproduce both ${{\omega }_{B}}t\ll 1$ [81–83] and ${{\omega }_{B}}t\gg 1$ limits [84, 125]; and addresses the general ${{\omega }_{B}}t\sim 1$ case, where formation of PSZS is expected [27] (see section 3).

5.3. Dyson equation for nearly periodic fluctuations

The same results that hold for equation (41) with broad and narrow frequency spectrum fluctuations can be obtained for equation (36), given the mapping of Equations (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⁷ . Furthermore, equation (36) can address the onset of nonlocal behaviors in EP transport, including the effects of equilibrium geometries and plasma nonuniformity discussed in section 2, which are expected to become progressively more important for increasing EP drive and nonperturbative response [27]. The present approach also shows that, even in the perturbative EP local limit and given the mapping of equations (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 ${{\bar{\omega }}_{d}}$ and ${{\hat{\omega }}_{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 the stochasticity threshold with all the necessary physics ingredients for a realistic comparison to 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 (see section 6).

Given a nearly periodic fluctuation $\delta {{\bar{\phi }}_{n}}(t)$ (spatial dependences are left implicit for the sake of notation simplicity), defined in equation (28) and with oscillation frequency centered about ${{\omega }_{0}}(\tau )={{\omega }_{0r}}(\tau )+{\rm i}{{\gamma }_{0}}(\tau )$ (slowly varying in time), we can write $\delta {{\bar{\phi }}_{n}}(t)\equiv {{{\rm lim} }_{\tau \to t}}\delta {{\bar{\varphi }}_{n}}(\tau ){\rm exp} \left( -{\rm i}{{\omega }_{0}}(\tau )t \right)$, with $\delta {{\bar{\varphi }}_{n}}(\tau )\equiv \delta {{\bar{\phi }}_{n0}}{\rm exp} \left[ -{\rm i}\int _{0}^{\tau }{{\omega }_{0}}(t^{\prime} ){\rm d}t^{\prime} +{\rm i}{{\omega }_{0}}(\tau )\tau \right]$. Thus, as shown in appendix C, $\delta {{\bar{\phi }}_{n}}(t)$ admits the Laplace transforms

$$\begin{eqnarray}&&\delta {{\widehat{{\bar{\phi }}}}_{n}}(\omega )=\frac{{\rm i}}{2\pi }\frac{\delta {{{\bar{\varphi }}}_{n}}(r,\tau )}{\omega -{{\omega }_{0}}(\tau )},\ \ \ \ \ {\rm and}\ \ \ \ \ \delta {{\widehat{{\bar{\phi }}}}_{-n}}(\omega )=\frac{{\rm i}}{2\pi }\frac{\delta {{{\bar{\varphi }}}_{-n}}(r,\tau )}{\omega +\omega _{0}^{*}(\tau )},\end{eqnarray} \tag{ 44 }$$

provided that $\mid t-\tau \mid \lesssim {{\tau }_{{\rm NL}}}\sim \mid {{\gamma }_{0}}{{\mid }^{-1}}$; and $\mid {{\dot{\omega }}_{0r}}\mid \ll \mid {{\gamma }_{0}}{{\omega }_{0r}}\mid$, $\mid {{\dot{\gamma }}_{0}}\mid \ll \mid \gamma _{0}^{2}\mid$. Therefore, equation (44) applies in general to adiabatic as well as nonadiabatic frequency chirping modes. When substituted back into equation (36), equation (44) yields the following form of the Dyson equation

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{F}}}_{0}}(\omega ) & = & \frac{{\rm i}}{\omega }{\rm St}{{{\hat{F}}}_{0}}(\omega )+\frac{{\rm i}}{\omega }{{{\hat{S}}}_{0}}(\omega )+\frac{{\rm i}}{2\pi \omega }{{{\bar{F}}}_{0}}(0)+\frac{e}{m}\frac{nc}{\omega ({\rm d}\psi /{\rm d}r)} \\ {} & {} & \times \frac{\partial }{\partial r}\left\{ \left[ \frac{Q_{n,{{\omega }_{0}}(\tau )}^{*}}{\omega _{0}^{*}(\tau )}\frac{{{{\hat{F}}}_{0}}\left( \omega -2{\rm i}{{\gamma }_{0}}(\tau ) \right)}{\omega -{{\omega }_{0}}(\tau )+n{{{\bar{\omega }}}_{d}}}+\frac{{{Q}_{n,{{\omega }_{0}}(\tau )}}}{{{\omega }_{0}}(\tau )} \right. \right. \\ {} & {} & \times \left. \left. \frac{{{{\hat{F}}}_{0}}\left( \omega -2{\rm i}{{\gamma }_{0}}(\tau ) \right)}{\omega +\omega _{0}^{*}(\tau )-n{{{\bar{\omega }}}_{d}}} \right]{{{\hat{\omega }}}_{dn}}\mid \delta {{{\bar{\varphi }}}_{n}}{{(r,\tau )}\mid ^{2}} \right\}, \\ \end{array}\end{eqnarray} \tag{ 45 }$$

where ${{Q}_{n,{{\omega }_{0}}(\tau )}}{{\hat{F}}_{\,0}}$ can be approximated by equation (37.)

Equation (45) is the renormalized form of PSZS generated by nearly periodic fluctuations, which allow investigating the self-consistent nonlinear evolution of SAW/DAW and EP PSZS when 'closed' by equations (18) (dropping the ZS shearing effect $\propto \;{{\partial }_{r}}{{\langle \delta {{L}_{g}}\rangle }_{z}}$) and (25). A crucial element, here, is the different argument of $\hat{F}$₀ on the lhs and rhs of equation (45), since it bears the information on the self-consistent interplay of nonlinear mode dynamics and EP transport. Assuming $\mid \omega \mid \gg \mid {{\gamma }_{0}}\mid \to 0$; i.e., assuming fixed amplitude fluctuations, allow us to solve equation (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 $\mid \omega \mid \sim \mid {{\gamma }_{0}}\mid$, and equation (45) describes much richer dynamics, including those to be discussed in section 6.

We analyze the nonlinear dynamics of a nearly periodic EPM burst using equations (18), (25) and (45). In particular, we consider a radially localized EPM wave packet with perpendicular wavelength ${{\lambda }_{\bot }}\ll {{L}_{{\rm 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 (see also section 3). Thus, ${{\theta }_{k0}}=0$ and $\partial D_{n}^{{\rm L}}/\partial {{\theta }_{k0}}=0$ in equation (25) [63, 64, 75, 76, 93]; and

$$\begin{eqnarray}&&\Lambda _{n}^{{\rm L}}\equiv {{\Lambda }_{{\rm T}}}=(1/2){{[({{\omega }^{2}}-\omega _{\ell }^{2})/(\omega _{u}^{2}-{{\omega }^{2}})]}^{1/2}}.\end{eqnarray} \tag{ 46 }$$

Here, the subscript T stands for 'toroidal'; and ${{\omega }_{\ell }}/{{\omega }_{u}}$ are the lower/upper SAW continuum accumulation point frequency. Using a simple $(s,\alpha )$ model for high aspect ratio (${{R}_{0}}/a\gg 1$) tokamak equilibria with circular shifted magnetic flux surfaces [131], with $s\equiv rq^{\prime} /q$ the magnetic shear and $\alpha \equiv -{{R}_{0}}{{q}^{2}}\beta ^{\prime}$ the 'ballooning' dimensionless pressure gradient parameter, it is possible to show [63, 64, 132]

$$\begin{eqnarray}&&\delta \bar{W}_{fn}^{L}=\frac{\mid s\mid \pi }{8}\left( 1+2\kappa (s)-\frac{\alpha }{{{\alpha }_{c}}} \right);\ \ \ \ \ {\rm and}\ \ \ \ \ \frac{1}{2}\frac{{{\partial }^{2}}D_{n}^{{\rm L}}}{\partial \theta _{k0}^{2}}=\frac{\mid s\mid \pi }{8}\kappa (s),\end{eqnarray} \tag{ 47 }$$

for $\mid s\mid ,\mid \alpha \mid \lt 1$ in equation (25), with ${{\alpha }_{c}}={{s}^{2}}/(1+\mid s\mid )$ and $\kappa (s)\simeq (1/2)(1+1/\mid s\mid ){\rm exp} (-1/\mid s\mid )$. Meanwhile, consistent with the 'fishbone' paradigm (see section 5), $\delta {{\bar{W}}_{kn}}$ can be written as [27, 63, 64]

$$\begin{eqnarray}\begin{array}{rcl} \delta {{{\bar{W}}}_{kn}} & = & \int \mathcal{E}{\rm d}\mathcal{E}{\rm d}\lambda \mathop{\sum }\limits_{{{v}_{\parallel }}/\mid {{v}_{\parallel }}\mid =\pm }\frac{{{\pi }^{2}}q{{R}_{0}}}{{{c}^{2}}k_{\vartheta }^{2}\mid s\mid }\frac{{{e}^{2}}}{m}\left( \frac{{{\tau }_{b}}{{n}^{2}}\bar{\omega }_{d}^{2}}{{{\omega }_{0}}(\tau )} \right) \\ {} & {} & \times \int _{-\infty }^{+\infty }\frac{\omega +{{\omega }_{0}}(\tau )}{n{{{\bar{\omega }}}_{d}}-{{\omega }_{0}}(\tau )-\omega }{{{\rm e}}^{-{\rm i}\omega t}}{{Q}_{n,{{\omega }_{0}}(\tau )}}{{{\hat{F}}}_{0}}(\omega ){\rm d}\omega , \\ \end{array}\end{eqnarray} \tag{ 48 }$$

where $\lambda =\mu {{B}_{0}}/\mathcal{E}$, ${{\tau }_{b}}=2\pi /{{\omega }_{b}}$, ${{k}_{\vartheta }}=-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 ${{\hat{\omega }}_{dn}}\simeq n{{\bar{\omega }}_{d}}$ in equation (35). The linear contribution $\delta \bar{W}_{kn}^{L}$ is readily obtained from equation (48) by substitution of ${{\hat{F}}_{\,0}}(\omega )={{(2\pi \omega )}^{-1}}{\rm i}{{\bar{F}}_{0}}(0)$, as suggested by equation (45). In this case, we assume an initial (equilibrium) EP isotropic slowing down distribution function of particles born at energy E_F. Thus,

$$\begin{eqnarray}&&{{\bar{F}}_{0}}(0)=\frac{3{{P}_{0{\rm EP}}}}{4\pi {{E}_{F}}}\frac{{\rm H}({{E}_{F}}/{{m}_{{\rm EP}}}-\mathcal{E})}{{{(2\mathcal{E})}^{3/2}}+{{(2{{E}_{c}}/{{m}_{{\rm EP}}})}^{3/2}}};\end{eqnarray} \tag{ 49 }$$

where ${\rm H}$ denotes the Heaviside step function, and the normalization condition is chosen such that the EP energy density is $(3/2){{P}_{0{\rm EP}}}$ for ${{E}_{F}}\gg {{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

$$\begin{eqnarray}&&\delta \bar{W}_{kn}^{L}({{\omega }_{0}})=\frac{3\pi {{(r/{{R}_{0}})}^{1/2}}}{8\sqrt{2}\mid s\mid }{{\alpha }_{{\rm EP}}}\left[ 1+\frac{{{\omega }_{0}}}{{{{\bar{\omega }}}_{dF}}}{\rm ln} \left( \frac{{{{\bar{\omega }}}_{dF}}}{{{\omega }_{0}}}-1 \right)+i\pi \frac{{{\omega }_{0}}}{{{{\bar{\omega }}}_{dF}}} \right],\end{eqnarray} \tag{ 50 }$$

where ${{\alpha }_{{\rm EP}}}=-8\pi {{R}_{0}}{{q}^{2}}P_{0{\rm EP}}^{\prime }/B_{0}^{2}$ and ${{\bar{\omega }}_{dF}}=n{{\bar{\omega }}_{d}}(\mathcal{E}={{E}_{F}}/{{m}_{{\rm EP}}})$. Considering a radially localized EP source, ${{\alpha }_{{\rm EP}}}={{\alpha }_{0{\rm EP}}}$ ${\rm exp} [-{{(r-{{r}_{0}})}^{2}}/L_{{\rm EP}}^{2}]$ [35]; and using equations (46), (47) and (50), it is possible to show that the linear EPM wave-packet, obtained as a solution of equation (25), has a characteristic radial envelope width $\sim \mid {{L}_{{\rm EP}}}/{{k}_{\vartheta }}{{\mid }^{1/2}}$, with ${{\lambda }_{\bot }}\sim \mid {{k}_{\vartheta }}{{\mid }^{-1}}\ll \mid {{L}_{{\rm EP}}}/{{k}_{\vartheta }}{{\mid }^{1/2}}\ll {{L}_{{\rm 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, $\propto {{\Lambda }_{T}}({{\omega }_{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 could adapt to the nonlinearly modified EP source; and, vice versa, that self-consistent evolution of EP, PSZS and EPM mode structures could take the form of avalanches [15, 35]. Thus, it is 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 the nonlinear dynamics of EP PSZS and EPM are dominated by $\delta \bar{W}_{kn}^{{\rm NL}}$. In fact, EP effects $\Lambda _{n}^{{\rm NL}}$ are negligible due to finite-orbit averaging in the kinetic/inertial layer, while $\delta \bar{W}_{fn}^{{\rm NL}}$ accounts for nonresonant EP responses [27, 63, 64]. Thus, in the present analysis, we can assume $\Lambda _{n}^{{\rm NL}}=\delta \bar{W}_{fn}^{{\rm NL}}=0$. We further set $S_{n}^{{\rm ext}}(r,t)=0$ in equation (25); and neglect source and collision terms in equation (45), as our aim is to investigate the nonlinear evolution of one single EPM burst excited by the initial unstable distribution given in equation (49). By direct substitution into equation (48) of ${{\hat{F}}_{0}}(\omega )$ from equation (45), and keeping the leading order terms in the time-scale ordering of equation (2), we obtain

$$\begin{eqnarray}\begin{array}{rcl} \delta \bar{W}_{kn}^{{\rm NL}} & \simeq & -\int \mathcal{E}{\rm d}\mathcal{E}{\rm d}\lambda \mathop{\sum }\limits_{{{v}_{\parallel }}/\mid {{v}_{\parallel }}\mid =\pm }\frac{{{\pi }^{2}}q{{R}_{0}}}{{{c}^{2}}k_{\vartheta }^{2}\mid s\mid }\frac{{{{\rm e}}^{2}}}{m}\int _{-\infty }^{+\infty }\frac{{{\tau }_{b}}{{n}^{2}}\bar{\omega }_{d}^{2}({{k}_{\vartheta }}/\Omega ){{{\rm e}}^{-{\rm i}\omega t}}}{(n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega )} \\ {} & {} & \times \frac{{{\partial }^{2}}}{\partial {{r}^{2}}}\left[ \frac{{{c}^{2}}k_{\vartheta }^{2}}{B_{0}^{2}}\frac{(2{\rm i}/\omega )({{\gamma }_{0}}+{\rm i}\omega )\mid \delta {{{\bar{\varphi }}}_{n}}{{(r,\tau )}\mid ^{2}}}{{{(n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})}^{2}}+{{({{\gamma }_{0}}+{\rm i}\omega )}^{2}}}\frac{\partial }{\partial r}{{{\hat{F}}}_{0}}(\omega -2{\rm i}{{\gamma }_{0}}) \right]{\rm d}\omega . \\ \end{array}\end{eqnarray} \tag{ 51 }$$

Note that, here, we have taken into account that nonlinear EPM dynamics occurs on meso-scales, intermediate between ${{\lambda }_{\bot }}\sim \mid {{k}_{\vartheta }}{{\mid }^{-1}}$ and L_EP, to have ${{\partial }_{r}}$ commute with quantities that vary on the equilibrium length scale.

In equation (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 $\propto 1/\omega$ dependence, while the other is due to the nearly singular/peaked structure of ${{\hat{F}}_{0}}(\omega -2{\rm i}{{\gamma }_{0}})$ at $\omega =2{\rm i}{{\gamma }_{0}}$. Due to the $\propto \;{\rm exp} (-{\rm i}\omega t)$ dependence of the Laplace integrand, the contribution from $\omega \simeq 0$ is exponentially smaller than that from $\omega \simeq 2{\rm i}{{\gamma }_{0}}$ in the growth and saturation phases of the burst; and thus is negligible. This further demonstrates the importance of the self-consistent interplay of nonlinear mode dynamics and EP transport, already noted in section 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 equation (51), at the leading order, we can rewrite

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial _{r}^{2}\left( \ldots \right)}{(n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega )} & = & \frac{{{\partial }^{2}}}{\partial {{r}^{2}}}\left[ \frac{\left( \ldots \right)}{(n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega )} \right] \\ {} & {} & +2\frac{(n{{\partial }_{r}}{{{\bar{\omega }}}_{d}}-{{{\dot{\omega }}}_{0r}}/{{v}_{{\rm g}}})-{\rm i}{{{\dot{\gamma }}}_{0}}/{{v}_{{\rm g}}}}{{{\left[ (n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega ) \right]}^{2}}}{{\partial }_{r}}\left( \ldots \right) \\ {} & {} & -2\frac{{{\left[ (n{{\partial }_{r}}{{{\bar{\omega }}}_{d}}-{{{\dot{\omega }}}_{0r}}/{{v}_{{\rm g}}})-{\rm i}{{{\dot{\gamma }}}_{0}}/{{v}_{{\rm g}}} \right]}^{2}}}{{{\left[ (n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega ) \right]}^{3}}}\left( \ldots \right) \\ {} & {} & +\frac{(n\partial _{r}^{2}{{{\bar{\omega }}}_{d}}-{{\ddot{\omega }}_{0r}}/v_{{\rm g}}^{2})-{\rm i}{{\ddot{\gamma }}_{0}}/v_{{\rm g}}^{2}}{{{\left[ (n{{{\bar{\omega }}}_{d}}-{{\omega }_{0r}})-{\rm i}({{\gamma }_{0}}-{\rm i}\omega ) \right]}^{2}}}\left( \ldots \right). \\ \end{array}\end{eqnarray} \tag{ 52 }$$

Noting that $\mid (n{{\bar{\omega }}_{d}}-{{\omega }_{0r}})\mid \sim \mid ({{\gamma }_{0}}-{\rm i}\omega )\mid \sim \mid {{k}_{r}}{{v}_{{\rm g}}}\mid \sim \tau _{{\rm NL}}^{-1}$; and consistent with the time scale ordering of equation (2) and discussed after equation (44), the dominant contributions on the rhs of equation (52) are given by the first three lines, where terms $\propto \;{{\dot{\gamma }}_{0}}$ and $\propto \;{{\ddot{\gamma }}_{0}}$ can be dropped. Meanwhile, given equation (A.35) for $\ell =0$, the contribution of second and third lines on the rhs are, respectively, $\sim \mid \tau _{{\rm NL}}^{2}{{\ddot{\Theta }}_{m,n,\ell =0}}\mid$ and $\sim \mid \tau _{{\rm NL}}^{4}{{\ddot{\Theta }}_{m,n,\ell =0}}{{\mid }^{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 (see section 5), account for wave-particle trapping [24, 25]. However, when phase locking occurs, equation (A.42), the second and third lines on the rhs of equation (52) become negligible, consistent with equation (6). Thus, phase locking de facto prevents wave-particle trapping from occurring (see appendix A.3). Note also that $\mid \tau _{{\rm NL}}^{2}{{\ddot{\Theta }}_{m,n,\ell =0}}\mid \sim {{\epsilon }_{{\dot{\omega }}}}\ll 1$ by definition, equation (A.45). Thus, on time scales shorter than $\epsilon _{{\dot{\omega }}}^{-1}{{\tau }_{{\rm NL}}}$, equation (52) with only the first line on the rhs can be used to calculate $\delta \bar{W}_{kn}^{{\rm NL}}$, which describes convective amplification of a nearly periodic phase-locked EPM burst when substituted back into Equation (25). Phase locking, equation (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 nonresonant EP responses as well as collisions and other relevant phenomena, neglected here but discussed in [15].

For the phase-locked EPM, the velocity space integral of equation (51) is dominated by resonant particles and by $\omega \simeq 2{\rm i}{{\gamma }_{0}}$, as discussed above. Thus, noting that $n{{\bar{\omega }}_{d}}\simeq -{{k}_{\vartheta }}\mathcal{E}/({{R}_{0}}\Omega )$ and ${{\tau }_{b}}\simeq 2\pi q{{R}_{0}}{{({{R}_{0}}/r)}^{1/2}}{{\mathcal{E}}^{-1/2}}$ for the simplified equilibrium considered in this section [64], equation (51) can be reduced to [27] ⁸

$$\begin{eqnarray}\begin{array}{rcl} \delta \bar{W}_{kn}^{{\rm NL}} & \simeq & -{\rm i}\frac{8{{\pi }^{4}}{{E}_{F}}}{\mid s\mid B_{0}^{2}}\frac{{{\omega }_{0r}}}{{{{\bar{\omega }}}_{dF}}}{{\left( \frac{r}{{{R}_{0}}} \right)}^{1/2}}{{q}^{2}}{{R}_{0}}\frac{{{\partial }^{2}}}{\partial {{r}^{2}}}\int _{-\infty }^{+\infty }\left[ \frac{{{c}^{2}}k_{\vartheta }^{2}}{B_{0}^{2}}\frac{\mid \delta {{{\bar{\phi }}}_{n}}(r,t){{\mid }^{2}}}{{{(2{{\gamma }_{0}}-{\rm i}\omega )}^{2}}} \right. \\ {} & {} & \times \left. \frac{\partial }{\partial r}{{\left( {{\mathcal{E}}^{3/2}}{{{\hat{F}}}_{0}}(\omega ) \right)}_{n{{{\bar{\omega }}}_{d}}={{\omega }_{0r}}}} \right]{{{\rm e}}^{-{\rm i}\omega t}}{\rm d}\omega =-{\rm i}\frac{8{{\pi }^{4}}{{E}_{F}}}{\mid s\mid B_{0}^{2}}\frac{{{\omega }_{0r}}}{{{{\bar{\omega }}}_{dF}}}{{\left( \frac{r}{{{R}_{0}}} \right)}^{1/2}} \\ {} & {} & \times {{q}^{2}}{{R}_{0}}\frac{{{\partial }^{2}}}{\partial {{r}^{2}}}\partial _{t}^{-2}\left[ \frac{{{c}^{2}}k_{\vartheta }^{2}}{B_{0}^{2}}\mid \delta {{{\bar{\phi }}}_{n}}{{(r,t)}\mid ^{2}}\frac{\partial }{\partial r}{{\left( {{\mathcal{E}}^{3/2}}{{F}_{0}}(t) \right)}_{n{{{\bar{\omega }}}_{d}}={{\omega }_{0r}}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 53 }$$

Furthermore, following [35] and noting that the analysis is limited to times shorter than $\epsilon _{{\dot{\omega }}}^{-1}{{\tau }_{{\rm NL}}}$, ${{F}_{\,0}}(t)\simeq {{\bar{F}}_{\,0}}(0)$ with ${{\bar{F}}_{\,0}}(0)$ given by equation (49). Thus, equation (51) can be further simplified and yields

$$\begin{eqnarray}&&\delta \bar{W}_{kn}^{{\rm NL}}\simeq {\rm i}\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]k_{\vartheta }^{2}\rho _{{\rm EP}}^{2}v_{{\rm EP}}^{2}\partial _{t}^{-2}\partial _{r}^{2}{{\mid \bar{A}}_{n0}}{{(r,t)}\mid ^{2}},\end{eqnarray} \tag{ 54 }$$

where $v_{{\rm EP}}^{2}={{E}_{F}}/{{m}_{{\rm EP}}}$, $\rho _{{\rm EP}}^{2}=v_{{\rm EP}}^{2}/\Omega _{{\rm EP}}^{2}$, $\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]$ is given by equation (50), and we have introduced the normalized amplitude [27]

$$\begin{eqnarray}&&({{e}_{{\rm EP}}}/{{E}_{F}})\delta {{\bar{\phi }}_{n}}(r,t)={{\bar{A}}_{n}}(r,t)={{{\rm e}}^{-{\rm i}{{\omega }_{0r}}t}}{{\bar{A}}_{n0}}(r,t).\end{eqnarray} \tag{ 55 }$$

Collecting terms, the desired nonlinear evolution equation for ${{\bar{A}}_{n0}}(r,t)$, equation (25), can then be rewritten as

$$\begin{eqnarray}\begin{array}{lll} {} & {} & \left[ {\rm i}\left( {{\Lambda }_{T}}({{\omega }_{0r}})-\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right] \right)-\left( \delta \bar{W}_{fn}^{{\rm L}}+\mathbb{R}{\rm e}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right] \right) \right]{{{\bar{A}}}_{n0}} \\ {} & {} & \quad =\left[ \frac{\partial \mathbb{R}{\rm e}\left[ \delta \bar{W}_{kn}^{L}({{\omega }_{0r}}) \right]}{\partial {{\omega }_{0r}}}\left( {\rm i}\frac{\partial }{\partial t} \right)+\frac{1}{2k_{\vartheta }^{2}\;{{s}^{2}}}\frac{{{\partial }^{2}}D_{n}^{L}}{\partial \theta _{k0}^{2}}\frac{{{\partial }^{2}}}{\partial {{r}^{2}}} \right. \\ {} & {} & \qquad \left. +{\rm i}\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]k_{\vartheta }^{2}\rho _{{\rm EP}}^{2}v_{{\rm EP}}^{2}\partial _{t}^{-2}\frac{{{\partial }^{2}}}{\partial {{r}^{2}}}{{{\mid \bar{A}}}_{n0}}\mid ^{2} \right]{{{\bar{A}}}_{n0}}, \\ \end{array}\end{eqnarray} \tag{ 56 }$$

where all contributions are defined in equations (46), (47) and (50).

Equation (56) readily recovers the linear limit [64, 93]. In the general case [27, 63], as suggested by the last term on the rhs, equation (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]

$$\begin{eqnarray}&&{{\bar{A}}_{n0}}(r,t)=U(\xi ){\rm exp} {{\int }^{t}}{{\gamma }_{0}}(t^{\prime} ){\rm d}t^{\prime} =W(\xi ){\rm exp} \left( {\rm i}\chi (\xi )+\int _{0}^{t}{{\gamma }_{0}}(t^{\prime} ){\rm d}t^{\prime} \right)\ ,\end{eqnarray} \tag{ 57 }$$

with

$$\begin{eqnarray}&&\xi \equiv {{k}_{n0}}\left( r-{{r}_{0}}-\int _{0}^{t}{{v}_{{\rm g}}}(t^{\prime} ){\rm d}t^{\prime} \right),\end{eqnarray} \tag{ 58 }$$

and k_n0 denoting the nonlinear radial wave vector. Noting that ${{\partial }_{r}}{{\bar{A}}_{n0}}={{k}_{n0}}{{\partial }_{\xi }}U(\xi ){\rm exp} {{\int }^{t}}{{\gamma }_{0}}(t^{\prime} ){\rm d}t^{\prime}$ and ${{\partial }_{t}}{{\bar{A}}_{n0}}=({{\gamma }_{0}}-{{k}_{n0}}{{v}_{{\rm g}}}{{\partial }_{\xi }})U(\xi ){\rm exp} {{\int }^{t}}{{\gamma }_{0}}(t^{\prime} ){\rm d}t^{\prime}$, there exist two regimes in the nonlinear EPM burst evolution: (i) ${{\gamma }_{0}}\gtrsim \mid {{k}_{n0}}{{v}_{{\rm g}}}\mid$, where the linear EPM is increasingly modified by the nonlinear interplay with EP transport as the fluctuation intensity increases; and (ii) ${{\gamma }_{0}}\lt \mid {{k}_{n0}}{{v}_{{\rm g}}}\mid$, where the nonlinear EPM-EP interaction dominates the dynamics. In this second case, it is possible to solve equation (56) in a simple analytical form. In fact, balancing the nonlinear term with the linear dispersiveness in equation (56) yields

$$\begin{eqnarray}&&{{v}_{{\rm g}}}={{\lambda }_{{\rm g}}}\hat{v}_{{\rm EP}}^{E\times B},\ \ \ \ \ {\rm and}\ \ \ \ \ k_{n0}^{2}=\frac{k_{\vartheta }^{2}}{\lambda _{{\rm g}}^{2}}\frac{{{s}^{2}}\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{L}({{\omega }_{0r}}) \right]}{{{\partial }^{2}}D_{n}^{L}/\partial \theta _{k0}^{2}},\end{eqnarray} \tag{ 59 }$$

with $\hat{v}_{{\rm EP}}^{E\times B}=(-{{k}_{\vartheta }}c/{{B}_{0}}){\rm max} [\delta {{\bar{\phi }}_{n}}(r,t)]$ the radial ${\boldsymbol{E}} \times {\boldsymbol{B}}$ velocity at the fluctuation peak, ${{\lambda }_{g}}$ an $\mathcal{O}(1)$ control parameter to be determined; and

$$\begin{eqnarray}&&\partial _{\xi }^{2}U={{\lambda }_{0}}U-2{\rm i}U\mid U{{\mid }^{2}},\end{eqnarray} \tag{ 60 }$$

where ${{\lambda }_{0}}\simeq -0.47+{\rm i}1.33$ corresponds to the ground state of the complex nonlinear oscillator described by equation (60), whose numerical solution is shown in figure 5, separating amplitude and phase, $U(\xi )=W(\xi ){{{\rm e}}^{{\rm i}\chi (\xi )}}$.

Zoom In Zoom Out Reset image size

With this solution, equation (56) becomes the 'nonlinear' EPM dispersion relation

$$\begin{eqnarray}&&D_{n}^{{\rm L}\ell }({{\omega }_{0}})-\frac{{{\lambda }_{0}}}{2\lambda _{{\rm g}}^{2}}\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]=0,\end{eqnarray} \tag{ 61 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} D_{n}^{{\rm L}\ell }({{\omega }_{0}}) & = & {\rm i}\left( {{\Lambda }_{T}}({{\omega }_{0r}})-\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right] \right)-\left( \delta \bar{W}_{fn}^{{\rm L}}+\mathbb{R}{\rm e}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right] \right) \\ {} & {} & -\left( {{\omega }_{0r}}\frac{\partial \mathbb{R}{\rm e}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]}{\partial {{\omega }_{0r}}} \right){\rm i}\frac{{{\gamma }_{0}}}{{{\omega }_{0r}}} \\ \end{array}\end{eqnarray} \tag{ 62 }$$

is the linearized 'local' EPM dispersion function [56], obtained neglecting the linear dispersiveness $\propto {{\partial }^{2}}D_{n}^{L}/\partial \theta _{k0}^{2}$ [64].

Equations (59) and (61) describe a one-parameter family, ${{\lambda }_{g}}$, of possible EPM wave-packets convectively amplified as they radially propagate with a group velocity $\propto \;\hat{v}_{{\rm EP}}^{E\times B}$. The dominant nonlinear mode is obtained for

$$\begin{eqnarray}&&\frac{{\rm d}{{\gamma }_{0}}}{{\rm d}\lambda _{{\rm g}}^{2}}=\frac{\partial {{\gamma }_{0}}}{\partial \lambda _{{\rm g}}^{2}}+\frac{\partial {{\gamma }_{0}}}{\partial {{\omega }_{0r}}}\frac{{\rm d}{{\omega }_{0r}}}{{\rm d}\lambda _{{\rm g}}^{2}}=0;\end{eqnarray} \tag{ 63 }$$

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 ${{\lambda }_{{\rm g}}}\simeq 0.5\div 0.6$ ⁹ . Thus, the EPM wave-packet propagates radially at a substantial fraction of the peak EP ${\boldsymbol{E}} \times {\boldsymbol{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 [16] for explaining EP secular loss due to 'fishbones'. Equation (61) applies at the instantaneous radial position ${{r}_{1}}(t)={{r}_{0}}+\int _{0}^{t}{{v}_{{\rm g}}}(t^{\prime} ){\rm d}t^{\prime}$ of the EPM wave-packet. The dependence on r₁, in the present model equilibrium, is explicit in terms $\propto {{\alpha }_{{\rm EP}}}$ and $\propto {{\bar{\omega }}_{dF}}$; and implicit in terms depending on ${{\omega }_{0}}$. Due to the $\propto {\rm ln} ({{\bar{\omega }}_{dF}}/{{\omega }_{0}}-1)$ term and ${{\omega }_{0}}/{{\bar{\omega }}_{dF}}$ dependences of $\delta \bar{W}_{kn}^{{\rm L}}$ in equation (50), the real EPM frequency is given by

$$\begin{eqnarray}&&{{\omega }_{0r}}/{{\bar{\omega }}_{dF}}({{r}_{1}})\simeq {\rm const},\end{eqnarray} \tag{ 64 }$$

at the leading oder, so that the phase-locking condition, equations (6) and (A.42), is readily satisfied, consistent with the ansatz and introductory discussions above in this section. Meanwhile, since $\mathbb{I}{\rm m}{{\lambda }_{0}}\gt 0$, the nonlinear EPM-EP interplay a strengthens the mode drive w.r.t. the linear value $\propto \mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]$. Thus, the convective amplification of the phase-locked EPM is an 'avalanche'[35]; and wave-packet amplification can continue until quenched by radial nonuniformities; i.e., the combined effect of decreasing values of linear drive $\mathbb{I}{\rm m}\left[ \delta \bar{W}_{kn}^{{\rm L}}({{\omega }_{0r}}) \right]\propto {{\alpha }_{{\rm EP}}}$ and increasing values of SAW continuum damping $\propto {{\Lambda }_{T}}({{\omega }_{0r}})$ (see section 2.2 and appendix A.3; and [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 nonadiabatic 'autoresonance' effect (see section 2.2), caused by the nonperturbative EP response and the ensuing kinetic EPM dispersiveness. Plasma instability and nonperturbative 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 nonadiabatic frequency chirping, the relevant nonlinear time scale is ${{\tau }_{{\rm NL}}}\sim \mid \partial _{t}^{-1}\mid \sim \mid {{A}_{n0}}{{\mid }^{-1}}$, consistent with equation (59). Furthermore, since ${{v}_{{\rm g}}}={{\lambda }_{{\rm g}}}\hat{v}_{{\rm EP}}^{E\times B}$, the EPM intensity is proportional to the square of the traveled distance, similar to the propagation of an FEL short optical pulse in the superradiant regime [62]. Similar to super-radiance 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 'super-radiance' in fusion plasmas [27, 134]. Other interesting analogies with neighboring research areas are further discussed in [27]. For example, equation (60) shows interesting analogies with $\partial _{\xi }^{2}U=U-2{{U}^{3}}$ yielding $U={\rm sech}(\xi )$; 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, equation (60) maintains its peculiarities, mostly due to its complex nature coming from the self-consistent interplay of EPM radial structure and EP transport.

Phase space zonal structures play fundamental roles in EP transport induced by SAW/DAW instabilities with $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \ll 1$ in burning-fusion plasmas. In this work, we have derived general equations that can be used as a 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 nonperturbative EP population in fusion plasmas with realistic equilibrium geometry and plasma nonuniformities. As a particular case of the general equations, we have derived a reduced description that could 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 'super-radiance' 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, $\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gt \mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}{{\mid }_{c}}$ ($\mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}{{\mid }_{c}}\lesssim {{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 nonlinear saturation dynamics, until it eventually dominates over resonance detuning. This transition is demonstrated in numerical simulation experiments of meso-scale 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 us to understand the onset of nonlocal behaviors, such as self-consistent interplay of mode structures and EP transport.

The strength of instability drive or, more precisely, the perturbative vs. non-perturbative 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 could evolve into hole/clump pairs with adiabatically sweeping frequency ($\mid \dot{\omega }\mid \ll \omega _{B}^{2}$), which are maintained at 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 nonperturbative, kinetic SAW/DAW dispersiveness, self-consistently modified by EP transport, causes nonadiabatic frequency evolution ($\mid \dot{\omega }\mid \sim \omega _{B}^{2}$). 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 nonperturbative EP population, is characterized by the peculiar system geometry and equilibrium non-uniformity. However, it also sheds light on physics aspects that could be of relevance to neighboring fields of research, such as space plasma physics, nonlinear dynamics, condensed matter and accelerator physics.

This work was supported by European Unions Horizon 2020 Research and Innovation program under grant agreement number 633053 as Enabling Research Project CfP-WP14-ER-01/ENEA_Frascati-01; and by US DoE, ITER-CN and NSFC No. 11235009 grants.

For gyrokinetic theory in axisymmetric toroidal systems considered in this work, two obvious pairs of action-angle coordinates for charged particle motions are $(m\mu ,\alpha )$, with $\mu =v_{\bot }^{2}/(2{{B}_{0}})+\ldots$ the magnetic moment (see, e.g., [60, 138]) and α the gyrophase; and $({{P}_{\phi }},\phi )$, with ${{P}_{\phi }}$ the canonical toroidal angular momentum and ϕ the toroidal angle¹⁰ . Here, perpendicular (⊥) and parallel (∥) directions are defined with respect to the unit vector ${\boldsymbol{b}} ={{{\boldsymbol{B}} }_{0}}/{{B}_{0}}$; and the axisymmetric equilibrium magnetic field configuration is assumed in the form

$$\begin{eqnarray}&&{{{\boldsymbol{B}} }_{0}}=F(\psi )\nabla \phi +\nabla \phi \times \nabla \psi ;\end{eqnarray} \tag{ A.1 }$$

$2\pi \psi$ denoting the poloidal magnetic flux within the magnetic surface labeled by ψ. Consistent with the main text, we adopt field-aligned toroidal flux coordinates $(r,\theta ,\zeta )$, with $r=r(\psi )$ a radial-like flux variable, the poloidal angle θ is chosen such that the Jacobian $\mathcal{J}={{(\nabla \psi \times \nabla \theta \cdot \nabla \zeta )}^{-1}}$ satisfies the condition $\mathcal{J}B_{0}^{2}$ being a magnetic flux function [139, 140]; and the general toroidal angle $\zeta \equiv \phi -\nu (\psi ,\theta )$, $\nu (\psi ,\theta )$ being a suitable periodic function of θ [141], such that the safety factor q,

$$\begin{eqnarray}&&q=\frac{{{{\boldsymbol{B}} }_{0}}\cdot \nabla \zeta }{{{{\boldsymbol{B}} }_{0}}\cdot \nabla \theta }=q(\psi ),\end{eqnarray} \tag{ A.2 }$$

is a function of ψ. At the leading order, ${{P}_{\phi }}$ is given by

$$\begin{eqnarray}&&{{P}_{\phi }}=\frac{e}{c}\left( F(\psi )\frac{{{v}_{\parallel }}}{\Omega }-\psi \right);\end{eqnarray} \tag{ A.3 }$$

with $\Omega =e{{B}_{0}}/(mc)$ the cyclotron frequency. The third pair of action-angle coordinates is $(J,{{\theta }_{{\rm c}}})$, with J the 'second invariant' and ${{\theta }_{{\rm c}}}$ the respective conjugate canonical angle

$$\begin{eqnarray}&&J=\oint {{v}_{\parallel }}{\rm d}l,\ \ \ \ \ {{\theta }_{c}}={{\omega }_{b}}\int _{0}^{\theta }{\rm d}\theta ^{\prime} /\dot{\theta }^{\prime} .\end{eqnarray} \tag{ A.4 }$$

Here, integrals are taken along the particle orbit in the equilibrium magnetic field, $dl$ denoting the arc-length element along it, and we have introduced the unified notation ${{\omega }_{b}}$ for bounce and transit frequency of trapped and circulating particles, respectively,

$$\begin{eqnarray}&&{{\omega }_{b}}(\mu ,J,{{P}_{\phi }})=\frac{2\pi }{\oint {\rm d}\theta /\dot{\theta }}.\end{eqnarray} \tag{ A.5 }$$

Given ${{H}_{0}}=m(v_{\parallel }^{2}+\mu {{B}_{0}})$ and the definitions of J and ${{\omega }_{b}}$,

$$\begin{eqnarray*}&&\frac{\partial J}{\partial {{H}_{0}}}=\frac{2\pi }{{{\omega }_{b}}}\ \ \leftrightarrow \ \ {{\omega }_{b}}=2\pi \frac{\partial {{H}_{0}}}{\partial J}.\end{eqnarray*}$$

Therefore, J can be readily computed given the reference magnetic equilibrium

$$\begin{eqnarray*}&&J=J(\mu ,{{P}_{\phi }},{{H}_{0}})\ \ \leftrightarrow \ \ {{H}_{0}}={{H}_{0}}(\mu ,J,{{P}_{\phi }}).\end{eqnarray*}$$

Equation (A.4) yields ${{\theta }_{{\rm c}}}={{\omega }_{b}}\tau$, with τ a time-like parameter tracking the particle position along the closed poloidal orbit of period $2\pi /{{\omega }_{b}}$. Thus, for a particle with given constants of motion,

$$\begin{eqnarray}&&r=\bar{r}+{{\tilde{\rho }}_{c}}({{\theta }_{c}}),\end{eqnarray} \tag{ A.6 }$$

where $\bar{r}$ and ${{\tilde{\rho }}_{{\rm c}}}({{\theta }_{{\rm c}}})$ are parameterized by $(\mu ,J,{{P}_{\phi }})$ and the initial particle radial position¹¹ ; and 'tilde' denotes a generic periodic function in ${{\theta }_{{\rm c}}}$ with zero average, which can be computed from the particle equations of motion in the equilibrium ${\boldsymbol{B}}$₀. Similarly, it is readily demonstrated that, for circulating particles (i.e.; particles whose ${{v}_{\parallel }}$ does not change sign)

$$\begin{eqnarray}&&\theta ={{\theta }_{{\rm c}}}+{{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}}),\end{eqnarray} \tag{ A.7 }$$

while, for magnetically trapped particles (i.e.; particles whose ${{v}_{\parallel }}$ changes sign twice along a closed orbit in the $(r,\theta )$ plane),

$$\begin{eqnarray}&&\theta ={{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}}).\end{eqnarray} \tag{ A.8 }$$

Note, again, that ${{\tilde{\Theta }}_{c}}$ is a periodic function of ${{\theta }_{c}}$ with zero average, which is parameterized by $(\mu ,J,{{P}_{\phi }})$. Thus, it takes different values in equations (A.7) and (A.8). At last, considering equations (A.2) and the definition of $\zeta \equiv \phi -\nu (\psi ,\theta )$, the parametric representation of ζ along the particle orbit in the equilibrium ${\boldsymbol{B}}$₀ is

$$\begin{eqnarray}&&\zeta ={{\bar{\omega }}_{d}}\tau +\bar{q}\theta +{{\tilde{\Xi }}_{{\rm c}}}({{\theta }_{{\rm c}}}).\end{eqnarray} \tag{ A.9 }$$

Here, θ is given by equations (A.7) or (A.8), ${{\bar{\omega }}_{d}}(\mu ,J,{{P}_{\phi }})$ is the toroidal precession frequency and

$$\begin{eqnarray}&&\bar{q}(\mu ,J,{{P}_{\phi }})=\frac{\oint q{\rm d}\theta }{\oint {\rm d}\theta }=\frac{\oint q\dot{\theta }{\rm d}{{\theta }_{{\rm c}}}}{\oint \dot{\theta }{\rm d}{{\theta }_{{\rm c}}}}.\end{eqnarray} \tag{ A.10 }$$

Note that, even though q does not depend on θ, it depends on the particle position because of the particle radial displacement. Meanwhile, from equation (A.9) and (A.10), one readily derives the following definition for the toroidal precession frequency

$$\begin{eqnarray}&&{{\bar{\omega }}_{{\rm d}}}(\mu ,J,{{P}_{\phi }})=\frac{{{\omega }_{b}}}{2\pi }\oint \left( \dot{\zeta }-q\dot{\theta } \right)\frac{{\rm d}\theta }{{\dot{\theta }}}.\end{eqnarray} \tag{ A.11 }$$

A.1. Wave-particle resonance condition

Using equations (A.6) to (A.9), the Fourier decomposition

$$\begin{eqnarray}&&f(r,\theta ,\zeta )=\mathop{\sum }\limits_{m,n\in \mathbb{Z}}{{{\rm e}}^{{\rm i}n\zeta -{\rm i}m\theta }}{{f}_{m,n}}(r)\end{eqnarray} \tag{ A.12 }$$

of a scalar function $f(r,\theta ,\zeta )$, describing a generic fluctuating field for which time dependence has been left implicit, can be projected along the unperturbed particle motion in the equilibrium ${\boldsymbol{B}}$₀ using the time-like parameter $\tau ={{\theta }_{c}}/{{\omega }_{b}}$; expressing $(r,\theta ,\zeta )$ in terms of the particle position along its unperturbed trajectory. More specifically, we have

$$\begin{eqnarray}&&f(r,\theta ,\zeta )=\mathop{\sum }\limits_{m,n,\ell \in \mathbb{Z}}{{{\rm e}}^{{\rm i}\left( n{{{\bar{\omega }}}_{d}}+\ell {{\omega }_{b}} \right)\tau }}{{\mathcal{P}}_{m,n,\ell }}\circ {{f}_{m,n}}(\bar{r}).\end{eqnarray} \tag{ A.13 }$$

Here, we have introduced the projection operators along the particle motion in the equilibrium ${\boldsymbol{B}}$₀, defined as

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{P}}_{m,n,\ell }}\circ {{f}_{m,n}}(\bar{r}) & \equiv & \frac{{{\lambda }_{m,n}}}{2\pi }\oint {\rm exp} \left\{ {\rm i}n{{\tilde{\Xi }}_{c}}({{\theta }_{c}})+{\rm i}\left[ n\bar{q}(\bar{r})-m \right] \right. \\ {} & {} & \times \left. {{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}}) \right\}{{f}_{m,n}}(\bar{r}+{{\tilde{\rho }}_{{\rm c}}}({{\theta }_{{\rm c}}})){{{\rm e}}^{-{\rm i}\ell {{\theta }_{{\rm c}}}}}d{{\theta }_{{\rm c}}}. \\ \end{array}\end{eqnarray} \tag{ A.14 }$$

Furthermore, ${{\lambda }_{m,n}}=1$ for magnetically trapped particles, while, for circulating particles,

$$\begin{eqnarray}&&{{\lambda }_{m,n}}={\rm exp} \left[ {\rm i}\left( n\bar{q}(\bar{r})-m \right){{\omega }_{b}}\tau \right].\end{eqnarray} \tag{ A.15 }$$

Note that equation (A.13) represents the value of $f(r,\theta ,\zeta )$ effectively experienced by the particle along its orbit, parameterized by the actions $(\mu ,J,{{P}_{\phi }})$, or equivalently by $(\mu ,{{P}_{\phi }},{{H}_{0}})$ (see above), and by the time-like variable $\tau ={{\theta }_{{\rm c}}}/{{\omega }_{b}}$. The integer $\ell \in \mathbb{Z}$ stands for the 'bounce harmonic', similarly to $m,n\in \mathbb{Z}$, denoting poloidal and toroidal Fourier harmonics. Equation (A.13) is completely defined, together with the functions ${{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}})$, ${{\tilde{\Xi }}_{{\rm c}}}({{\theta }_{{\rm c}}})$ and ${{\tilde{\rho }}_{{\rm c}}}({{\theta }_{{\rm c}}})$ that are obtained from the integration of particle motion, once the reference equilibrium ${\boldsymbol{B}}$₀ 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, equation (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 is thus 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,\theta ,\zeta )$ are left implicit to simplify notation. Once a monochromatic wave $f(r,\theta ,\zeta )\propto {\rm exp} (-{\rm i}{{\omega }_{0}}t)$ is assumed, the resonance condition, i.e., the stationarity of wave-particle phase in the particle-moving frame, where $\tau \leftrightarrow t$, is readily derived from equation (A.13) and yields

$$\begin{eqnarray}&&{{\omega }_{0}}={{\omega }_{0}}(\mu ,J,{{P}_{\phi }})=n{{\bar{\omega }}_{d}}+\ell {{\omega }_{b}}\end{eqnarray} \tag{ A.16 }$$

for magnetically trapped particles, while for circulating particles

$$\begin{eqnarray}&&{{\omega }_{0}}={{\omega }_{0}}(\mu ,J,{{P}_{\phi }})=n{{\bar{\omega }}_{d}}+\ell {{\omega }_{b}}+\left( n\bar{q}(\bar{r})-m \right){{\omega }_{b}}.\end{eqnarray} \tag{ A.17 }$$

Here, we remind readers that ${{\omega }_{b}}$ stands for both bounce and transit frequencies. We also note that ${{\bar{\omega }}_{d}}$ is typically negligible for circulating particles, except for particles that are close to the trapped-to-circulating boundary in the action space.

A.2. The effect of fluctuations

When considering the effect of fluctuations, for every bounce/transit; i.e., a $2\pi$ increment in the angle ${{\theta }_{{\rm c}}}$, a shift will be accumulated in the wave-particle phase ('resonance detuning'), because of radial nonuniformities and the dependence of ${{\bar{\omega }}_{d}}$ and ${{\omega }_{b}}$ on $J,{{P}_{\phi }}$ that are no longer conserved in the nonlinear regime (μ is conserved for low-frequency fluctuations, analyzed here, for which $\mid {{\omega }_{0}}\mid \ll \mid \Omega \mid$). Moreover, the wave-particle phase in the mode structure decomposition in action-angle variables will also be shifted, since ${{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}})$, ${{\tilde{\Xi }}_{{\rm c}}}({{\theta }_{{\rm c}}})$ and ${{\tilde{\rho }}_{{\rm c}}}({{\theta }_{{\rm c}}})$ are not simple periodic functions of ${{\theta }_{{\rm 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\gg \mid {{\omega }_{0}}{{\tau }_{{\rm NL}}}{{\mid }^{-1}}\sim \mid {{\gamma }_{{\rm L}}}/{{\omega }_{0}}\mid \gg {{\epsilon }_{\omega }}\equiv \mid {{\omega }_{0}}/\Omega \mid$ (see section 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 equation (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 $\sim {{\epsilon }_{\delta }}\equiv \mid \delta {{{\boldsymbol{B}} }_{\bot }}/{{B}_{0}}\mid \ll 1$ [60], yielding the expressions for $\delta \dot{r}$ ($\delta \dot{\psi }$), $\delta \dot{\theta }$ and $\delta \dot{\zeta }$, as well as for $\dot{H}$₀, ${{\dot{P}}_{\phi }}$ and $\dot{J}$. Here, $(\delta r,\delta \theta ,\delta \zeta )$ denotes the change in particle position due to the small but finite amplitude fluctuations, which are also the cause of $\delta {{H}_{0}}$, $\delta {{P}_{\phi }}$ and $\delta J$ and finite $\dot{H}$ a₀, ${{\dot{P}}_{\phi }}$ and $\dot{J}$. 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\left\langle \delta {{\phi }_{g}}-{{v}_{\parallel }}\delta {{A}_{\parallel g}}/c \right\rangle -\bar{H}$ (see, e.g., [70, 142]), whose conservation can be readily cast as

$$\begin{eqnarray}&&\frac{\partial }{\partial t}{{\dot{\Pi }}_{\phi }}+\frac{\partial }{\partial \zeta }\mathop{{\bar{H}}}\limits^{.}=0.\end{eqnarray} \tag{ A.18 }$$

Here, $\bar{H}$ is the the gyrocenter Hamiltonian up to order $\sim {{\epsilon }_{\delta }}$ [60]

$$\begin{eqnarray}&&\bar{H}={{H}_{0}}+e\left\langle \delta {{\phi }_{g}} \right\rangle -e\frac{{{v}_{\parallel }}}{c}\left\langle \delta {{A}_{\parallel g}} \right\rangle ;\end{eqnarray} \tag{ A.19 }$$

and ${{\Pi }_{\phi }}$ is the extension of the canonical toroidal angular momentum, equation (A.3), in the presence of fluctuations and for ${\boldsymbol{B}}$₀ expressed as in equation (A.1),

$$\begin{eqnarray}&&{{\Pi }_{\phi }}={{P}_{\phi }}+\frac{e}{c}\frac{F(\psi )}{{{B}_{0}}}\left\langle \delta {{A}_{\parallel g}} \right\rangle .\end{eqnarray} \tag{ A.20 }$$

Furthermore, it is readily shown that, for trapped particles, we have

$$\begin{eqnarray}&&\theta ={{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}})+\Delta \theta ,\end{eqnarray} \tag{ A.21 }$$

while, for circulating particles,

$$\begin{eqnarray}\begin{array}{rcl} \theta & = & {{\omega }_{b}}\tau +{{\tilde{\Theta }}_{{\rm c}}}({{\theta }_{{\rm c}}})+\Delta \theta \\ {} & {} & +\frac{\partial {{\omega }_{b}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} +\frac{\partial {{\omega }_{b}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} , \\ \end{array}\end{eqnarray} \tag{ A.22 }$$

with

$$\begin{eqnarray}&&\Delta \theta =\int _{0}^{\tau }\delta \dot{\theta }\frac{{\rm d}\theta }{{\dot{\theta }}}.\end{eqnarray} \tag{ A.23 }$$

Note that the difference between equations (A.21) and (A.22) stems from the different nature of unperturbed orbits, respectively described by equations (A.8) and (A.7), and from the varied expression of ${{\theta }_{{\rm c}}}\to {{\theta }_{{\rm c}}}+\delta {{\theta }_{{\rm c}}}$ ($\delta {{\theta }_{{\rm c}}}=\int _{0}^{\tau }\delta {{\omega }_{b}}{\rm d}\tau ^{\prime}$) for circulating particles in the presence of fluctuations. Similarly, the radial motion is described as

$$\begin{eqnarray}&&r=\bar{r}+{{\tilde{\rho }}_{{\rm c}}}({{\theta }_{{\rm c}}})+\Delta r,\end{eqnarray} \tag{ A.24 }$$

$$\begin{eqnarray}&&\Delta r=\int _{0}^{\tau }\delta \dot{r}\frac{{\rm d}\theta }{{\dot{\theta }}},\end{eqnarray} \tag{ A.25 }$$

while ζ is given by

$$\begin{eqnarray}\begin{array}{rcl} \zeta & = & {{{\bar{\omega }}}_{d}}\tau +\bar{q}\theta +{{\tilde{\Xi }}_{{\rm c}}}({{\theta }_{{\rm c}}})+\Delta \zeta \\ {} & {} & +\frac{\partial {{{\bar{\omega }}}_{d}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} +\frac{\partial {{{\bar{\omega }}}_{d}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} \\ {} & {} & +\bar{q}\left( \frac{\partial {{\omega }_{b}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} +\frac{\partial {{\omega }_{b}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} \right)+{{\omega }_{b}}\frac{{\rm d}\bar{q}}{{\rm d}\bar{r}}\int _{0}^{\tau }\delta r{\rm d}\tau ^{\prime} , \\ \end{array}\end{eqnarray} \tag{ A.26 }$$

$$\begin{eqnarray}&&\Delta \zeta =\int _{0}^{\tau }\delta \dot{\zeta }\frac{{\rm d}\theta }{{\dot{\theta }}}\ .\end{eqnarray} \tag{ A.27 }$$

Here, the terms involving ${{\omega }_{b}}$ and its derivatives with respect to ${{P}_{\phi }}$ and J are to be considered for circulating particles only. In fact, these terms are originated by the effect of fluctuations on the term $\propto \bar{q}{{\theta }_{c}}=\int _{0}^{\tau }{{\omega }_{b}}\bar{q}{\rm d}\tau ^{\prime}$, due to the unperturbed circulating particle motion, equation (A.7), which has no correspondence in the unperturbed trapped particle motion, equation (A.8). In particular, terms $\propto \partial {{\omega }_{b}}/\partial {{P}_{\phi }},\partial {{\omega }_{b}}/\partial J$ have the same origin of those included in equation (A.22) for circulating particles and excluded in equation (A.21) for magnetically trapped particles. Similarly, the term involving $\sim {{\omega }_{b}}({\rm d}\bar{q}/{\rm d}\bar{r})$ stems from $\int _{0}^{\tau }{{\omega }_{b}}\delta \bar{q}{\rm d}\tau ^{\prime}$ for circulating particles only.

With the above results, we can extend the mode structure decomposition in action-angle variables, equation (A.13), to include nonlinear particle orbit distortions due to fluctuations:

$$\begin{eqnarray}&&f(r,\theta ,\zeta )=\mathop{\sum }\limits_{m,n,\ell \in \mathbb{Z}}\lambda _{m,n}^{{\rm NL}}{{{\rm e}}^{{\rm i}\left( n{{{\bar{\omega }}}_{d}}+\ell {{\omega }_{b}} \right)\tau +{\rm i}\Theta _{m,n,\ell }^{{\rm NL}}}}{{\mathcal{P}}_{m,n,\ell }}\circ {{{\rm e}}^{\Delta r{{\partial }_{{\bar{r}}}}}}{{f}_{m,n}}(\bar{r}),\end{eqnarray} \tag{ A.28 }$$

where ${{\mathcal{P}}_{m,n,\ell }}\mid {{{\rm e}}^{\Delta r{{\partial }_{{\bar{r}}}}}}{{f}_{m,n}}$ is defined as in equation (A.14), $\lambda _{m,n}^{{\rm NL}}$ is the nonlinear extension of equation (A.15), i.e., $\lambda _{m,n}^{{\rm NL}}=1$ for trapped particles, while for circulating particles

$$\begin{eqnarray}\begin{array}{rcl} \lambda _{m,n}^{{\rm NL}} & = & {\rm exp} \left[ {\rm i}\left( n\bar{q}(\bar{r})-m \right)\left( \frac{\partial {{\omega }_{b}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} \right. \right. \\ {} & {} & \left. \left. +\frac{\partial {{\omega }_{b}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} \right)+{\rm i}n{{\omega }_{b}}\frac{{\rm d}\bar{q}}{{\rm d}\bar{r}}\int _{0}^{\tau }\delta r{\rm d}\tau ^{\prime} \right] \\ \end{array}\end{eqnarray} \tag{ A.29 }$$

and the definition of $\Theta _{m,n,\ell }^{{\rm NL}}$ is

$$\begin{eqnarray}\begin{array}{rcl} \Theta _{m,n,\ell }^{{\rm NL}} & = & n\Delta \zeta -m\Delta \theta +n\left( \frac{\partial {{{\bar{\omega }}}_{{\rm d}}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} +\frac{\partial {{{\bar{\omega }}}_{d}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} \right) \\ {} & {} & +\ell \left( \frac{\partial {{\omega }_{b}}}{\partial {{P}_{\phi }}}\int _{0}^{\tau }\delta {{P}_{\phi }}{\rm d}\tau ^{\prime} +\frac{\partial {{\omega }_{b}}}{\partial J}\int _{0}^{\tau }\delta J{\rm d}\tau ^{\prime} \right). \\ \end{array}\end{eqnarray} \tag{ A.30 }$$

Equation (A.28) can be rewritten in a more compact form, introducing the definition

$$\begin{eqnarray}&&{{\Theta }_{m,n,\ell }}\equiv \Theta _{m,n,\ell }^{{\rm NL}}-{\rm i}{\rm ln} \lambda _{m,n}^{{\rm NL}}-\int _{0}^{\tau }\Delta \omega (t){\rm d}t,\end{eqnarray} \tag{ A.31 }$$

so that the nonlinear frequency shift $\Delta \omega (t)=\omega (t)-{{\omega }_{0}}$ for a nearly monochromatic wave (see section 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 waves still applies for frequency chirping modes as long as $\mid \dot{\omega }\mid =\mid \Delta \dot{\omega }\mid \ll \mid {{\gamma }_{{\rm L}}}{{\omega }_{0}}\mid$; i.e., for both adiabatic ($\mid \dot{\omega }\mid \ll \omega _{B}^{2}$) and nonadiabatic ($\mid \dot{\omega }\mid \lesssim \omega _{B}^{2}$) frequency sweeping, ${{\omega }_{B}}$ denoting the wave-particle trapping frequency. Meanwhile, we have, at the lowest order,

$$\begin{eqnarray}&&{{\mathcal{P}}_{m,n,\ell }}\circ {{{\rm e}}^{\Delta r{{\partial }_{r}}}}{{f}_{m,n}}(\bar{r})\simeq {{\mathcal{P}}_{m,n,\ell }}\mid {{f}_{m,n}}(\bar{r}+\Delta r),\end{eqnarray} \tag{ A.32 }$$

with $\Delta r$ being the bounce-averaged nonlinear radial particle displacement, accumulated over many bounces/transits ($2\pi$ increments in ${{\theta }_{c}}$). For each bounce/transit, the bounce-averaging involved in the definition of ${{\mathcal{P}}_{m,n,\ell }}\mid {{f}_{m,n}}(\bar{r})$, equation (A.14), is carried out at the 'effective' $\bar{r}$ location shifted by the accumulated bounce-averaged nonlinear radial particle displacement. Therefore, equation (A.28) can be finally cast as

$$\begin{eqnarray}&&f(r,\theta ,\zeta )=\mathop{\sum }\limits_{m,n,\ell \in \mathbb{Z}}{{{\rm e}}^{{\rm i}\left( n{{{\bar{\omega }}}_{d}}+\ell {{\omega }_{b}} \right)\tau +{\rm i}{{\Theta }_{m,n,\ell }}}}{{\mathcal{P}}_{m,n,\ell }}\circ {{f}_{m,n}}(\bar{r}+\Delta r).\end{eqnarray} \tag{ A.33 }$$

Thus, the dominant effect of nonlinear physics on wave-particle resonances enters via the nonlinear phase shift ${{\Theta }_{m,n,\ell }}$ ('resonance detuning') and/or the ${{\mathcal{P}}_{m,n,\ell }}\circ {{f}_{m,n}}(\bar{r}+\Delta 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' versus 'radial decoupling' in numerical simulations of nonlinear dynamics of SAW excited by EPs, is given in [53].

A.3. Finite interaction length and finite interaction time

By direct inspection of equations (A.29) and (A.30), 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 equation (A.18), one generally has $\Delta r/r\sim \Delta {{P}_{\phi }}/{{P}_{\phi }}\sim ({{\omega }_{*{\rm EP}}}/{{\omega }_{0}})\Delta {{H}_{0}}/{{H}_{0}}$ [27, 106], with ${{\omega }_{*{\rm EP}}}$ the EP diamagnetic frequency. Thus, not only radial decoupling, but also resonance detuning for EPs, is expected to be dominated by nonlinear radial displacement, since $\mid {{\omega }_{*E}}/{{\omega }_{0}}\mid \gg 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., $\ell =\mathcal{O}(1)$. Specializing, for simplicity, to the case of a toroidal equilibrium with major radius R₀ and shifted circular magnetic surfaces, we readily obtain

$$\begin{eqnarray}&&{{\dot{\Theta }}_{m,n,l}}\simeq n\frac{{\rm d}\bar{q}}{{\rm d}\bar{r}}{{\omega }_{t}}\Delta r-\Delta \omega ,\end{eqnarray} \tag{ A.34 }$$

for circulating particles, having denoted the transit frequency as ${{\omega }_{t}}$ and neglected ${{\bar{\omega }}_{d}}$; whereas, for magnetically trapped particle resonance,

$$\begin{eqnarray}&&{{\dot{\Theta }}_{m,n,l}}\simeq \left( n\frac{\partial {{{\bar{\omega }}}_{d}}}{\partial \bar{r}}+\ell \frac{\partial {{{\bar{\omega }}}_{b}}}{\partial \bar{r}} \right)\Delta r-\Delta \omega .\end{eqnarray} \tag{ A.35 }$$

Meanwhile, for both circulating and magnetically trapped particles, the equation for $\Delta r$ is

$$\begin{eqnarray}\begin{array}{rcl} \Delta \dot{r} & \simeq & {\rm i}\frac{c}{{{B}_{0}}}\mathop{\sum }\limits_{m,n,\ell \in \mathbb{Z}}\frac{n\bar{q}}{{\bar{r}}}{{{\rm e}}^{{\rm i}\left( n{{{\bar{\omega }}}_{d}}+\ell {{\omega }_{b}} \right)\tau +{\rm i}{{\Theta }_{m,n,\ell }}}} \\ {} & {} & \times {{\mathcal{P}}_{m,n,\ell }}\circ {{\left\langle \frac{\bar{r}q}{r\bar{q}}\left[ 1-\frac{{{v}_{\parallel }}}{q{{R}_{0}}\omega }(nq-m) \right]\delta {{\phi }_{g}} \right\rangle }_{m,n}}(\bar{r}+\Delta r), \\ \end{array}\end{eqnarray} \tag{ A.36 }$$

having used $\delta {{A}_{\parallel m,n}}\simeq (nq-m)/(q{{R}_{0}})(c/\omega )\delta {{\phi }_{m,n}}$; i.e., $\delta {{E}_{\parallel }}\simeq 0$ at the lowest order for SAW. For a resonance given by equations (A.16) or (A.17), equations (A.34) to (A.36) could describe an isolated resonance as well as resonance overlap in the case of many modes. It is also worthwhile to emphasize that, for the same toroidal mode number n, different poloidal mode numbers m are generally coherent, since they belong to the same (generally nonlinear) mode structure. From equations (A.34) to (A.36) with $\Delta \omega =0$, the wave-particle trapping frequency ${{\omega }_{B}}$ is readily estimated as

$$\begin{eqnarray}&&\omega _{B}^{2}\simeq {{\lambda }_{n}}\mid \frac{\omega }{r}\delta {{\mathop{{\boldsymbol{\bar{X}}} }\limits^{.}}_{\bot }}\mid \simeq {{\lambda }_{n}}\mid \frac{\omega }{r}\frac{nq}{r}\frac{c}{{{B}_{0}}}\delta \phi \mid ,\end{eqnarray} \tag{ A.37 }$$

with ${{\lambda }_{n}}=\mid nrq^{\prime} \mid$ for circulating, ${{\lambda }_{n}}=1$ for magnetically trapped particles; and $\delta {{\mathop{{\boldsymbol{\bar{X}}} }\limits^{.}}_{\bot }}$ the fluctuating ${\boldsymbol{E}} \times {\boldsymbol{B}}$ guiding-center drift velocity [60, 138]. In the derivation of equation (A.37), the estimate ${{d}_{tt}}\Delta r\sim \omega _{B}^{2}\Delta r$ has been adopted, and resonant particles have been assumed; i.e., Equation (A.16) for magnetically trapped and equation (A.17) for circulating particles, respectively. Given ${{\omega }_{B}}$, the separatrix width in particle phase space, $\Delta {{r}_{sx}},$ is estimated as $\Delta {{r}_{sx}}\simeq \mid \delta {{\mathop{{\boldsymbol{\bar{X}}} }\limits^{.}}_{\bot }}/{{\omega }_{B}}\mid$. Thus, $\Delta {{r}_{sx}}$ and ${{\omega }_{B}}$ are of order $\sim \epsilon _{\delta }^{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 nonuniform system. For $\Delta \omega =0$; i.e., no frequency chirping, the finite interaction time for resonant particles defines the characteristic nonlinear time scale and is given by

$$\begin{eqnarray}&&{{\dot{\Theta }}_{m,n,\ell }}{{\tau }_{{\rm NL}}}\sim 1;\ \ \ \ \ {\rm and}\ \ \ \ \ {{\tau }_{{\rm NL}}}\sim 1/(3{{\gamma }_{{\rm L}}}).\end{eqnarray} \tag{ A.38 }$$

Note that the $\gamma _{{\rm L}}^{-1}$ enters via the anticipated saturation process (either resonance detuning or radial decoupling, or both, with varying relative weight); i.e., saturation occurs when ${{\tau }_{{\rm NL}}}\sim 1/(3{{\gamma }_{{\rm 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 AEs in toroidal plasmas [145]. Meanwhile, neglecting the radial mode structure, the finite interaction length, $\Delta {{r}_{L}}$, is obtained from equation (A.34) as

$$\begin{eqnarray}&&3{{\gamma }_{{\rm L}}}\sim nq^{\prime} \omega \Delta {{r}_{L}}\ \ \end{eqnarray} \tag{ A.39 }$$

for circulating particles; while, from equation (A.35),

$$\begin{eqnarray}&&3{{\gamma }_{{\rm L}}}\sim \omega (\Delta {{r}_{L}}/r)\end{eqnarray} \tag{ A.40 }$$

for magnetically trapped particles, where characteristic radial scale lengths of ${{\omega }_{t}}$, ${{\bar{\omega }}_{d}}$ and ${{\omega }_{b}}$ are taken to be $\sim r$. From these estimates, we readily recognize that the finite interaction length, $\Delta {{r}_{L}}$, for circulating resonant particles is typically smaller than the characteristic width, $\sim \mid nq^{\prime} {{\mid }^{-1}}$ [63, 64, 111], of one poloidal harmonic, $\delta {{\phi }_{m,n}}$; whereas it might be comparable or larger than $\sim \mid nq^{\prime} {{\mid }^{-1}}$ for magnetically trapped resonant particles if $({{\gamma }_{{\rm L}}}/\omega )\gtrsim \mid 3nrq^{\prime} {{\mid }^{-1}}$. For high-n modes as expected in ITER [63, 64], this condition is $({{\gamma }_{{\rm L}}}/\omega )\gtrsim {{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 are intrinsically nonlocal; i.e., characterized by meso-scales, larger than $\mid nq^{\prime} {{\mid }^{-1}}$ and shorter than the equilibrium scale length, and convective (ballistic) processes [27]. Different trapped particle versus circulating particle behaviors are very general, as they are connected with the properties of equations (A.34) and (A.35), and are expected to play an 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 [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 nonadiabatic [$\mid \dot{\omega }\mid \lesssim \omega _{B}^{2}$; cf. discussions following equation (A.31)] frequency sweeping as, e.g., in NSTX experiments with neutral beam injection [150–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 equations (A.34) to (A.36), which have the same structure of the original equations derived in [16]; and predict secular (ballistic) particle motions if the ${{\dot{\Theta }}_{m,n,\ell }}=0$ resonance condition is preserved as the particles are transported out of the system. In fact, equations (A.39) and (A.40) give 'mode particle pumping' by 'phase locking' for circulating and magnetically trapped particles when, respectively [16, 153],

$$\begin{eqnarray}&&\Delta \dot{\omega }\simeq n\frac{{\rm d}\bar{q}}{{\rm d}\bar{r}}{{\omega }_{t}}\Delta \dot{r},\end{eqnarray} \tag{ A.41 }$$

and

$$\begin{eqnarray}&&\Delta \dot{\omega }\simeq \left( n\frac{\partial {{{\bar{\omega }}}_{d}}}{\partial \bar{r}}+\ell \frac{\partial {{{\bar{\omega }}}_{b}}}{\partial \bar{r}} \right)\Delta \dot{r}.\end{eqnarray} \tag{ A.42 }$$

It is worthwhile noting that, because of equations (A.41) and (A.42) together with equations (A.36) and (A.37), phase locking and secular (ballistic) particle transport always imply nonadiabatic frequency sweeping, $\Delta \dot{\omega }\sim \omega _{B}^{2}.$ Another important feature of the phase-locking condition and $\Delta \dot{r}$ given by equation (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., [3, 151]), and in numerical simulations of, e.g., nonlinear EPM evolutions [73, 94–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 Equations (A.41) and (A.42) in the $\Delta \omega =0$ no-frequency sweeping limit, to, respectively,

$$\begin{eqnarray}&&3{{\gamma }_{{\rm L}}}\sim nq^{\prime} \omega {{\epsilon }_{{\dot{\omega }}}}\Delta {{r}_{{\rm L}}},\end{eqnarray} \tag{ A.43 }$$

and

$$\begin{eqnarray}&&3{{\gamma }_{{\rm L}}}\sim \omega {{\epsilon }_{{\dot{\omega }}}}(\Delta {{r}_{L}}/r).\end{eqnarray} \tag{ A.44 }$$

Here, ${{\epsilon }_{{\dot{\omega }}}}$ is defined as

$$\begin{eqnarray}&&{{\epsilon }_{{\dot{\omega }}}}\equiv \mid \frac{{{{\dot{\Theta }}}_{m,n,\ell }}(\Delta \omega )}{{{{\dot{\Theta }}}_{m,n,\ell }}(\Delta \omega =0)}\mid ,\end{eqnarray} \tag{ A.45 }$$

where numerator and denominator are computed from equation (A.31), respectively, at the actual value of the nonlinear frequency shift (frequency chirping) and at $\Delta \omega =0$ (no chirping). The value of ${{\epsilon }_{{\dot{\omega }}}}$ denotes the effectiveness of phase locking; it depends on the wave-dispersive properties [27] and, in general, ${{\epsilon }_{{\dot{\omega }}}}\lt 1$ requires the EP dynamics nonperturbative. A similar argument applies to the finite interaction time, which is also extended by $\sim \epsilon _{{\dot{\omega }}}^{-1}$. Because of the extended interaction length, circulating resonant particles can also be displaced by a significant fraction of $\sim \mid nq^{\prime} {{\mid }^{-1}}$ before resonance detuning sets in. Thus, in the presence of significant nonadiabatic 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–104].

The secular resonant-particle motion, predicted by equations (A.41) and (A.42) is restricted to a group of particles that have a 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 lose the resonance condition after the finite interaction time/length, they can be substituted by others; and the process can continue until quenched by equilibrium nonuniformity (see sections 2.2 and 6; and [15, 26, 27]). In fact, nonadiabatic frequency sweeping prevents the de facto wave-particle trapping to occur, and the particle can enter and leave the resonant region during the nonlinear evolution, which maximizes wave-particle power exchange as well as mode growth [27] (see section 6).

Here, expressions for ${{\Lambda }_{n}}$, $\delta {{\bar{W}}_{fn}}$ and $\delta {{\bar{W}}_{kn}}$ to be used in equation (23) are given, without derivation, for the readers' convenience. Detailed derivations of these expressions can be found in [27, 63, 64].

The generalized 'inertia', ${{\Lambda }_{n}}$, accounting for the radially local structures of the SAW continuous spectrum, is given by

$$\begin{eqnarray}&&{\rm i}{{\Lambda }_{n}}\equiv \frac{1}{2}{{\left( \delta \hat{\Psi }_{-n{{0}^{+}}}^{\dagger }\delta {{{\hat{\Psi }}}_{n{{0}^{+}}}} \right)}^{-1}}\left[ \delta \hat{\Psi }_{-n}^{\dagger }(\vartheta ){{\partial }_{\vartheta }}\delta {{{\hat{\Psi }}}_{n}}(\vartheta ) \right]_{\vartheta \to {{0}^{-}}}^{\vartheta \to {{0}^{+}}},\end{eqnarray} \tag{ B.1 }$$

where $\delta {{\hat{\Psi }}_{n{{0}^{+}}}}={{{\rm lim} }_{\vartheta \to {{0}^{+}}}}\delta {{\hat{\Psi }}_{n}}(\vartheta )$, and $\delta {{\hat{\Psi }}_{n}}(\vartheta )$ is obtained from the solution of the short radial scale limit ($k_{\bot }^{2}\simeq k_{r}^{2}$) of the gyrokinetic vorticity equation written as [27, 63, 64]

$$\begin{eqnarray}&&\begin{array}{l} \left( \frac{{{\partial }^{2}}}{\partial {{\vartheta }^{2}}}-\frac{\partial _{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}}{{{{\hat{\kappa }}}_{\bot }}} \right)\delta {{{\hat{\Psi }}}_{n}}-\frac{{{\mathcal{J}}^{2}}B_{0}^{2}}{v_{A}^{2}}\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial t}+{\rm i}{{\omega }_{*pi}} \right. \\ \quad -\left. \frac{3}{4}k_{\vartheta }^{2}\rho _{{\rm i}}^{2}\hat{\kappa }_{\bot }^{2}\left( \frac{\partial }{\partial t}+{\rm i}{{\omega }_{*pi}}+{\rm i}{{\omega }_{*Ti}} \right) \right]\delta {{{\hat{\Phi }}}_{n}} \\ \quad -\mathop{{\rm lim} }\limits_{{{{\hat{\kappa }}}_{\bot }}\to \infty }\frac{4\pi {{\mathcal{J}}^{2}}{{B}_{0}}}{ck_{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}}{\boldsymbol{b}} \times {\boldsymbol{ \kappa }} \cdot \nabla \sum {{\left\langle m\left( \mu {{B}_{0}}+v_{\parallel }^{2} \right){{J}_{0}}\frac{\partial }{\partial t}\delta {{{\hat{g}}}_{n}} \right\rangle }_{v}} \\ \quad +\left[ {\rm NL}\;{\rm TERMS} \right]=0. \\ \end{array}\end{eqnarray} \tag{ B.2 }$$

Here, we have adopted equation (21) for representing spatiotemporal structures of SAW/DAW, ${{{\boldsymbol{k}} }_{\bot }}=-{\rm i}{{\nabla }_{\bot }}$, $k_{\bot }^{2}\equiv \hat{\kappa }_{\bot }^{2}k_{\vartheta }^{2}$, ${{k}_{\vartheta }}\equiv -nq/r$, $\delta {{\hat{\Psi }}_{n}}={{\hat{\kappa }}_{\bot }}\delta {{\hat{\psi }}_{n}}$, $\delta {{\hat{\Phi }}_{n}}={{\hat{\kappa }}_{\bot }}\delta {{\hat{\phi }}_{n}}$, ${{\rho }_{i}}={{({{T}_{i}}/{{m}_{i}})}^{1/2}}/{{\Omega }_{i}}$ is the thermal ion Larmor radius, $\mathcal{J}$ is the Jacobian [cf. appendix A, following equation (A.1)], summation is extended to all particle species, and we have introduced the thermal ion diamagnetic frequencies

$$\begin{eqnarray}\begin{array}{rcl} {{\omega }_{*pi}} & = & {{\omega }_{*ni}}+{{\omega }_{*Ti}}, \\ {{\omega }_{*ni}} & = & {{\left( \frac{{{T}_{0}}c}{e{{n}_{0}}{{B}_{0}}} \right)}_{i}}({\boldsymbol{b}} \times \nabla {{n}_{0i}})\cdot {{{\boldsymbol{k}} }_{\bot }}, \\ {{\omega }_{*Ti}} & = & {{\left( \frac{c}{e{{B}_{0}}} \right)}_{i}}({\boldsymbol{b}} \times \nabla {{T}_{0i}})\cdot {{{\boldsymbol{k}} }_{\bot }}. \\ \end{array}\end{eqnarray} \tag{ B.3 }$$

Furthermore, for the sake of simplicity, we have assumed the long wavelength limit $k_{\bot }^{2}\rho _{i}^{2}\ll 1$, except in the magnetic curvature coupling term, where J₀ is maintained for brevity. Interested readers can find a discussion of the more general ${{k}_{\bot }}{{\rho }_{i}}\sim 1$ case in [27, 63, 64]. Note that the $\sim \delta {{\hat{g}}_{n}}$ term generally includes nonlinear particle response via $\delta {{\hat{g}}_{n}}$ itself. Finally, formally nonlinear terms, denoted as $\left[ {\rm NL}\;{\rm TERMS} \right]$, are given by

$$\begin{eqnarray}\begin{array}{rcl} \left[ {\rm NL}\;{\rm TERMS} \right] & = & \mathop{{\rm lim} }\limits_{{{{\hat{\kappa }}}_{\bot }}\to \infty }\frac{{{\mathcal{J}}^{2}}B_{0}^{2}}{k_{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}v_{{\rm A}}^{2}}\frac{{{\partial }^{2}}}{\partial {{t}^{2}}}{{\mathcal{P}}_{Bn}}\left( \frac{\delta {{\varrho }_{{\rm m}}}}{{{\varrho }_{{\rm m}0}}}\nabla _{\bot }^{2}\delta \phi \right), \\ {} & {} & -\mathop{{\rm lim} }\limits_{{{{\hat{\kappa }}}_{\bot }}\to \infty }\frac{{{\mathcal{J}}^{2}}{{B}_{0}}}{ck_{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}}{{\mathcal{P}}_{Bn}}\left( {\boldsymbol{b}} \times \nabla \delta {{A}_{\parallel }}\cdot \nabla \nabla _{\bot }^{2}\delta {{A}_{\parallel }} \right) \\ {} & {} & -\mathop{{\rm lim} }\limits_{{{{\hat{\kappa }}}_{\bot }}\to \infty }\mathop{\sum }\limits_{\ne e}\frac{2\pi e{{\mathcal{J}}^{2}}B_{0}^{2}}{ck_{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}{{\Omega }^{2}}}{{\mathcal{P}}_{Bn}}\left[ {\boldsymbol{b}} \times \nabla \left( \nabla _{\bot }^{2}\delta \phi \right)\cdot \nabla {{\left\langle \mu \delta g \right\rangle }_{v}} \right. \\ {} & {} & -\left. {\boldsymbol{b}} \times \nabla \delta \phi \cdot \nabla {{\left\langle \mu \nabla _{\bot }^{2}\delta g \right\rangle }_{v}}-\nabla _{\bot }^{2}\left( {\boldsymbol{b}} \times \nabla \delta \phi \cdot \nabla {{\left\langle \mu \delta g \right\rangle }_{v}} \right) \right]. \\ \end{array}\end{eqnarray} \tag{ B.4 }$$

Here, $\delta {{\varrho }_{{\rm m}}}$ stands for plasma mass density fluctuation about the equilibrium value ϱ_m0, summation is extended to all particle species except thermal electrons and, in order to simplify notation, we have introduced the operators ${{\mathcal{P}}_{Bn}}$ [63, 64],

$$\begin{eqnarray}&&{{\mathcal{P}}_{Bn}}(r,\vartheta )\;:\;f(r,\theta ,\zeta )\mapsto {{\hat{f}}_{n}}(r,\vartheta ),\end{eqnarray} \tag{ B.5 }$$

which act on a generic function $f(r,\theta ,\zeta )$ in $(r,\theta ,\zeta )$ space, producing the corresponding ${{\hat{f}}_{n}}(r,\vartheta )$ component in $(r,\vartheta )$ space, consistent with equation (21). In equation (B.4), the first line represents the polarization current nonlinearity [130], the second one accounts for the Maxwell stress, while the third line is the kinetic expression representing both nonlinear diamagnetic response and Reynolds stress for $k_{\bot }^{2}\rho _{i}^{2}\ll 1$ [27, 63].

The expressions of $\delta {{\bar{W}}_{fn}}$ and $\delta {{\bar{W}}_{kn}}$ are obtained by construction of a quadratic form starting from the gyrokinetic vorticity equation [27, 63], and are given by

$$\begin{eqnarray}\begin{array}{rcl} \delta {{{\bar{W}}}_{nf}}+\delta {{{\bar{W}}}_{nk}} & = & {{\left( \delta \hat{\Phi }_{-n{{0}^{+}}}^{\dagger }\delta {{{\hat{\Phi }}}_{n{{0}^{+}}}} \right)}^{-1}}\frac{1}{2}\int _{-\infty }^{\infty }\left[ {{\left( \frac{\partial }{\partial \vartheta }\delta {{{\hat{\Phi }}}_{-n}} \right)}^{\dagger }}\left( \frac{\partial }{\partial \vartheta }\delta {{{\hat{\Phi }}}_{n}} \right) \right. \\ {} & {} & +\frac{\partial _{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}}{{{{\hat{\kappa }}}_{\bot }}}\delta \hat{\Phi }_{-n}^{\dagger }\delta {{{\hat{\Phi }}}_{n}}+\delta \hat{\Phi }_{-n}^{\dagger }\frac{{{\mathcal{J}}^{2}}B_{0}^{2}}{v_{{\rm A}}^{2}}\frac{\partial }{\partial t}\left( \frac{\partial }{\partial t}+{\rm i}{{\omega }_{*pi}} \right)\delta {{{\hat{\Phi }}}_{n}} \\ {} & {} & +\left. \delta \hat{\Phi }_{-n}^{\dagger }\frac{4\pi {{\mathcal{J}}^{2}}{{B}_{0}}}{ck_{\vartheta }^{2}{{{\hat{\kappa }}}_{\bot }}}{\boldsymbol{b}} \times {\boldsymbol{ \kappa }} \cdot \nabla \sum {{\left\langle m\left( \mu {{B}_{0}}+v_{\parallel }^{2} \right){{J}_{0}}\frac{\partial }{\partial t}\delta {{{\hat{g}}}_{n}} \right\rangle }_{v}} \right]{\rm d}\vartheta . \\ \end{array}\end{eqnarray} \tag{ B.6 }$$

Here, we noted that $\delta {{\hat{\Phi }}_{n}}(\vartheta )\simeq \delta {{\hat{\Psi }}_{n}}(\vartheta )$ for regular SAW/DAW radial structures away from resonances with the SAW continuous spectrum. Furthermore, the nonlinear contribution to $\delta {{\bar{W}}_{nf}}+\delta {{\bar{W}}_{nk}}$ is mainly due to the $\sim \delta {{\hat{g}}_{n}}$ term, while formally nonlinear terms in the gyrokinetic vorticity equation [cf. Equation (B.4)] predominantly contribute to ${{\Lambda }_{n}}$ [27, 63]. The nonlinearity due to EP response in the $\sim \delta {{\hat{g}}_{n}}$ term is what is considered in this work for the analyses of EPM avalanches, carried out in section 6 adopting the 'fishbone' paradigm (see section 5). The expressions of $\delta {{\bar{W}}_{fn}}$ and $\delta {{\bar{W}}_{kn}}$ in equation (B.6) can be extracted, dividing $\delta {{\hat{g}}_{n}}$ in 'fluid' and 'kinetic' responses as specified in equation (31) [17]. The expression for $\delta {{\bar{W}}_{kn}}$ is then obtained from equation (B.6) using the 'kinetic' response only in the $\sim \delta {{\hat{g}}_{n}}$ term. Note that the remaining $\delta {{\bar{W}}_{fn}}$ does not include the (linear) 'kink drive' contribution, which is typically negligible for SAW/DAW and can be anyway readily included for application of the present theoretical framework to MHD modes [27, 63, 64].

In this Appendix, we justify the Laplace representation of equation (44) for nearly periodic fluctuations, and discuss its validity in terms of relevant time scales.

We introduce the slow and fast time variables τ and t, with $\tau \simeq t$ but $\mid {{\partial }_{t}}\mid \gg \mid {{\partial }_{\tau }}\mid$ and ${{d}_{t}}={{\partial }_{t}}+{{\partial }_{\tau }}$. Furthermore, given a nearly periodic fluctuation $\delta {{\bar{\phi }}_{n}}(t)$ (spatial dependences are left implicit for the sake of notation simplicity), defined in equation (28) and with oscillation frequency centered about ${{\omega }_{0}}(\tau )={{\omega }_{0r}}(\tau )+{\rm i}{{\gamma }_{0}}(\tau )$ (slowly varying in time), we let

$$\begin{eqnarray}\begin{array}{rcl} \delta {{{\bar{\phi }}}_{n}}(t) & \equiv & \mathop{{\rm lim} }\limits_{\tau \to t}\delta {{\varphi }_{n}}(t,\tau )\equiv \mathop{{\rm lim} }\limits_{\tau \to t}\delta {{{\bar{\varphi }}}_{n}}(\tau ){\rm exp} \left( -{\rm i}{{\omega }_{0}}(\tau )t \right) \\ \delta {{{\bar{\varphi }}}_{n}}(\tau ) & \equiv & \delta {{{\bar{\phi }}}_{n0}}{\rm exp} \left[ -{\rm i}\int _{0}^{\tau }{{\omega }_{0}}(t^{\prime} ){\rm d}t^{\prime} +{\rm i}{{\omega }_{0}}(\tau )\tau \right]. \\ \end{array}\end{eqnarray} \tag{ C.1 }$$

Noting that

$$\begin{eqnarray*}&&({{\partial }_{t}}+{{\partial }_{\tau }})\delta {{\varphi }_{n}}(t,\tau )=\left[ -{\rm i}{{\omega }_{0}}(\tau )+{\rm i}{{{\dot{\omega }}}_{0}}(\tau )(\tau -t) \right]\delta {{\varphi }_{n}}(t,\tau ),\end{eqnarray*}$$

equation (C.1) admits the Laplace transforms

$$\begin{eqnarray}&&\delta {{\widehat{{\bar{\phi }}}}_{n}}(\omega )=\frac{{\rm i}}{2\pi }\frac{\delta {{{\bar{\varphi }}}_{n}}(r,\tau )}{\omega -{{\omega }_{0}}(\tau )},\ \ \ \ \ {\rm and}\ \ \ \ \ \delta {{\widehat{{\bar{\phi }}}}_{-n}}(\omega )=\frac{{\rm i}}{2\pi }\frac{\delta {{{\bar{\varphi }}}_{-n}}(r,\tau )}{\omega +\omega _{0}^{*}(\tau )},\end{eqnarray} \tag{ C.2 }$$

provided that $\mid t-\tau \mid \lesssim {{\tau }_{{\rm NL}}}\sim \mid {{\gamma }_{0}}{{\mid }^{-1}}$; and $\mid {{\dot{\omega }}_{0r}}\mid \ll \mid {{\gamma }_{0}}{{\omega }_{0r}}\mid$, $\mid {{\dot{\gamma }}_{0}}\mid \ll \mid \gamma _{0}^{2}\mid$. Thus, considering the time-scale ordering of equation (2), equation (C.2) applies in general to adiabatic as well as nonadiabatic frequency chirping modes; and is adopted in this work as equation (44).