Evidence of"Two Plasmon"decay of energetic particle induced geodesic acoustic mode

Secondary low frequency mode generation by energetic particle induced geodesic acoustic mode (EGAM) observed in LHD experiment is studied using nonlinear gyrokinetic theory. It is found that the EGAM frequency can be significantly higher than local geodesic acoustic mode (GAM) frequency in low collisionality plasmas, and it can decay into two GAMs as its frequency approaches twice GAM frequency, in a process analogous to the well-known two plasmon decay instability. The condition for this process to occur is also discussed.


Introduction
Confinement of plasmas in magnetically confined fusion devices [1] is one of the key issues for the sustained burning and fusion gain. The anomalous transport induced by micro-scale turbulence excited by expansion free energy intrinsic to confinement is thus an important topic of fusion plasma research [2,3]. Zonal field structures (ZFS) with toroidally and nearly poloidally symmetric mode structures (n = 0/m ≃ 0, ±1 · · ·, with n/m being the toroidal/poloidal mode numbers, respectively) are generally recognized to regulate micro-scale drift wave turbulence (DW) including drift Alfvén waves (DAWs) and their associated transport by scattering into stable short radial wavelength regimes [4][5][6][7][8]. Interested readers may refer to a recent work [9] discussing the role of ZFS and phase space zonal structures (PSZS) [10,11] as the generator of nonlinear equilibria with (suppressed) turbulence [6,12].
Geodesic acoustic modes (GAMs) [13,14] are ZFS unique in toroidal plasmas; oscillating at a finite frequency in the ion sound frequency regime due to toroidicity induced plasma compression. GAMs are predominantly electrostatic radial corrugations with the scalar potential characterized by n = 0/m ≃ 0, ±1, · · ·, and an up-down anti-symmetric m ≃ 1 density perturbation. Consequently, GAM, being toroidally symmetric, cannot by itself induce cross-field transport. Instead, GAMs can regulate DW/DAWs and the associated transport via spontaneous excitation, since they can scatter DW/DAWs into stable short radial wavelength domain [14][15][16][17]. Interested readers may refer to Ref. [18] for a brief review of kinetic theories of GAM, including linear dispersion relation, excitation by super-thermal energetic particles (EPs) and nonlinear interaction with DWs/DAWs. Due to their finite frequencies, GAMs can resonantly interact and exchange energy with charged particles including super-thermal EPs, leading to, respectively, collisionless Landau damping [19,20] and resonant excitation by EPs [21][22][23]. The excitation of EP-induced GAM (EGAM) was first analytically investigated in Ref. [24], showing the dominant role played by wave-particle resonant interactions and free energy associated with EP velocity space anisotropy, with application to different scenarios characterized by different EP distribution functions [25][26][27][28][29][30]. The global features due to GAM continuum coupling are also investigated [31], yielding a finite threshold due to continuum damping, or mode conversion to propagating kinetic GAM, as kinetic dispersiveness of both thermal and energetic ions are properly accounted for [26,32]. The EGAMs are typically characterized by a global mode structure with frequency lower than local GAM frequency and radial extension determined by EP density profile and kinetic dispersiveness, L ∼ √ L n ρ h , with L n being the EP density profile scale length, and ρ h being EP characteristic orbit width. Nonlinear saturation of EGAM due to wave-particle trapping [33] in the weak drive limit is investigated [34,35], and the EGAM induced EP loss via pitch angle scattering into unconfined orbits is presented in Ref. [36].
Recently, a peculiar phenomena was reported in Large Helical Device (LHD) low density experiments with neutral beam injection (NBI) heating [23]. During the discharge with relatively high plasma temperature (central electron temperature up to ∼ 7keV) and low plasma density (line averaged density ∼ 0.1 × 10 19 m −3 ), the EGAM frequency is observed to be significantly higher than local GAM frequency. The interpretation was given in Ref. [30], where it was shown that, due to the low collisionality in the high-temperature/low-density discharges, the injected beam ions are not fully slowed down, and form a bump-on-tail type EP distribution function. As a consequence, the free energy associated with the positive slope in the low energy side of the distribution function, provides an additional drive, that generates a new unstable branch with the frequency being significantly higher than local GAM frequency [30]. A similar interpretation was also provided in Ref. [29] considering a similar distribution function induced by low charge exchange rate.
It was further found [37] that, as the high frequency EGAM (will be denoted as "primary EGAM" in the rest of the paper for apparent reasons) chirped up to twice local GAM frequency, possibly as a result of the EGAM induced pitch angle scattering to lower pitch angle and thus higher v domain, a "secondary mode" could be strongly driven unstable, with the frequency close to the local GAM frequency [37]. The experimental observations are nicely recovered in a MHD-kinetic hybrid simulation using MEGA code [38,39], and it is found that, besides the secondary mode strong excitation, the EPs "driving" the secondary mode are the same as those linearly driving the primary EGAM, though the secondary mode frequency is only half of that of the primary mode. It is also found that the secondary mode generation still persists when the "fluid nonlinearity" is turned off; and, as the primary mode frequency keeps chirping up, the secondary mode frequency is almost unchanged, suggesting it is a normal mode of the system itself. One speculation by Ref. [39] is that the secondary mode is driven by "nonlinear resonance" [40,41], which, however, is typically associated with finite amplitude fluctuations, and is not satisfied in the condition for the onset of the secondary mode.
In this work, we will show that, the mechanism underlying the secondary mode generation observed in Ref. [37] is analogous to the well-known two plasmon decay process [42,43], where an incident electromagnetic wave decays into two plasma waves identical to each other. The primary mode corresponds to the linearly unstable high frequency EGAM given in Ref. [30], while the secondary mode corresponds to a linearly stable branch described by the same linear dispersion relation, with the frequency determined by the local GAM frequency. The secondary mode can be nonlinearly driven unstable as the primary mode frequency approaches twice the secondary mode frequency (which is close to the local GAM frequency), which minimizes the threshold due to frequency mismatch. The interpretation is consistent with the crucial elements from experimental observation [37] and numerical simulation [39]; thus, it provides and illuminates the underlying physics picture.
The rest of the paper is organized as follows. In Sec. 2, the theoretical model is given. The linear particle response to EGAM is derived in Sec. 3, where the linear EGAM properties in the LHD low collisionality plasma are also reviewed. In Sec. 4, the nonlinear decay of the primary mode into two low frequency GAMs is investigated, and its correspondence to the experimental observations and MEGA simulations are discussed. Finally, a brief summary and discussion is presented in Sec. 5. For the self-containedness of the materials, the properties of the three branches of the linear dispersion relation are briefly reproduced in Appendix A, while the frequency upchirping of the primary EGAM, being an important condition for the decay process, is briefly outlined in Appendix B.

Theoretical Model
In the above described nonlinear decay process, a primary EGAM (Ω 0 ≡ Ω 0 (ω 0 , k r,0 )) decays into two secondary modes with almost identical frequency, Ω 1 ≡ Ω 1 (ω 1 , k r,1 ) and Ω 2 ≡ Ω 2 (ω 2 , k r,2 ). All the three modes, are EGAM/GAMs with n = 0/m ≃ 0, thus the toroidal/poloidal wavenumber matching condition is naturally satisfied, and only the constraint on frequency and radial wavenumber matching condition is needed. For the modes with predominantly electrostatic polarization, the equations describing the nonlinear decay of the linearly unstable primary EGAM, can be derived from the charge quasi-neutrality condition. Assuming for simplicity the bulk ions and EPs are the same ion species with unit electric charge, the quasi-neutrality condition can be written as: with the perturbed particle density derived from the perturbed distribution function as Here, subscripts s = e, i, h denote electrons, ions and energetic (hot) particles, respectively, F 0,s is the equilibrium distribution function, ∂ E F 0,s is the short notation for ∂F 0,s /∂E, with E ≡ v 2 /2, J 0 (k ⊥ ρ L,s ) is the Bessel function of zero index accounting for finite Larmor radius (FLR) effects, k ⊥ is the perpendicular wavenumber, ρ L ≡ v ⊥ /Ω s is the Larmor radius with Ω s being gyro-frequency, and · · · denotes velocity space integration. δH is the non-adiabatic particle response to GAM/EGAM, and can be derived from nonlinear gyrokinetic equation [44]: Here, ω tr ≡ v /(qR 0 ) is the transit frequency, ω d = k r v d is the magnetic drift frequency associated with geodesic curvature, and v d = (v 2 ⊥ + 2v 2 ) sin θ/(2ΩR 0 ) ≡ v d sin θ is the magnetic drift velocity, J k ≡ J 0 (k ⊥ ρ L ) for simplicity of notation, and k=k ′ +k ′′ denotes the usual selection rule for frequency/wavenumber matching condition for the nonlinear mode coupling. The other notations are standard.
The thermal plasma linear response to GAM/EGAM, can be readily obtained from Ref. [20], and after surface averaging, equation (1) with the leading order linear thermal plasma response properly accounted for by the first term, from which the linear dispersion relation of GAM can be obtained, and (· · ·) ≡ 2π 0 (· · ·)dθ/(2π) is the surface averaging. Furthermore, ω G ≡ 7/4 + τ v it /R is the leading order GAM frequency with τ = T e /T i being the temperature ratio and v it being ion thermal velocity, while higher order terms such as FLR and/or finite parallel compression can be straightforwardly accounted for by replacing ω G with more accurate expression [20]. The linear EGAM dispersion relation can be derived by keeping δn L h in equation (4), with the free energy from EP velocity space anisotropy and the characteristic features of the EGAM determined by the specific EP distribution function [24-26, 29, 30]. Nonlinear modulation of GAM/EGAM by thermal plasma nonlinearity, can be accounted for by δn N L i , including excitation by DW/DAWs [14,16,45]. Here, super-scripts "L" and "NL" represent linear and nonlinear responses, respectively.
For the LHD low collisionality discharge [23] of interest, where plasma density is very low while the electron temperature is very high, the EP is characterized by a not fully slowed down distribution function, analogous to a bump-on-tail case, and the corresponding linear properties of the EGAM are investigated in Ref. [30]. The linear EGAM properties are the basis of the present nonlinear analysis, and, thus, for the self-containedness of the present nonlinear analysis, the results of Ref. [30] will be briefly summarized in Sec. 3 In this work, both the thermal plasma and EP induced nonlinear coupling are consistently derived [46,47], by including their nonlinear contribution to density perturbation in the quasi-neutrality condition. As we will show a posteriori, the nonlinear coupling is dominated by the EP finite orbit width effects (FoWs), with resonant EPs playing the crucial role [47,48]. The thermal plasma contribution to the nonlinear coupling [46], meanwhile, will be shown to be negligible compared to the dominant role of EPs.

Linear properties
For the completeness of this work, we will briefly derive the linear EP response to GAM/EGAM, which will be used in deriving the nonlinear response of EP to the secondary EGAMs. Separating the linear from nonlinear responses by taking δH k = δH L k + δH N L k , the linear EP response to EGAM can be derived from the linear gyrokinetic equation, Equation (5) can be solved by transforming into the EP drift orbit center coordinate by taking δH L k ≡ e iΛ k δH L dk , with Λ k satisfying ω tr ∂ θ Λ k +ω d,k = 0. Here, for simplicity of discussion and focus on proof of principle demonstration, well circulating EPs are assumed, thus, variation of v along the magnetic field is neglected, and one has Λ k =Λ k cos θ, withΛ k = k rρd being radial wave-number normalized to drift orbit width, andρ d =v dr /ω tr is the EP magnetic drift orbit width. The generalization to finite inverse aspect ratio case as well as general geometry and particle orbits is straightforward. Furthermore, τ ≡ T e /T i ≪ 1 is assumed for simplicity, such that one has δφ G ≃ δφ G while ω tr,e ≫ ω G is still satisfied ‡. We then have, from which δH L dk can be derived as Here, l is integer, and l ≡ ∞ l=−∞ will be used later for simplicity of notation. In deriving equation (7), the Jacobi-Anger expansion e iΛ k cos θ = l i l J l (Λ k )e ilθ is used. Pulling back to the EP guiding center coordinate, we then have, the linear well-circulating EP response to EGAM Interested readers may refer to Ref. [20] for the contribution of finite τ to linear GAM dispersion relation via the m = 0 components of the scalar potential.
Substituting equation (8) into the surface averaged quasi-neutrality condition, equation (4), one then obtains, e m i n 0 k 2 with the EGAM dispersion relation given by Here, only the l = ±1 transit harmonic are kept in equation (10), in consistency with the typicalΛ k ≪ 1 ordering for EGAMs with global mode structure. The EGAM stability depends sensitively on the specific EP equilibrium distribution function. For the not fully slowed down EPs in LHD experiments [23] due to low collisionality, the distribution can be modelled as which can be solved as a dynamic solution of the Fokker-Planck equation with the collisional operator dominated by thermal electron induced slowing down [30]. Here, is the time dependent lower end of the distribution function, and E crit is the critical energy at which the EP pitch angle scattering rate off thermal ions is comparable with the slowing down rate γ c [49]. The normalization to EP density n b gives ). Note that, c 0 is proportional to Γ/γ c with Γ being the NBI particle flux and γ c is the slowing-down rate on thermal electrons. The resulting dispersion relation is given as [30] Here, is the ratio of EP to bulk plasma density, are the transit frequencies defined at E b and E L , respectively. The first term in the square bracket corresponds to the slowing-down distribution induced logarithmic singularity, and it is destabilizing for λ 0 B 0 > 2/5 [26]; while the single pole like singularity in the second term is from the low energy side cutoff, and it is always destabilizing [30].
The dispersion relation (12) have three branches, i.e., an unstable branch determined by and close to ω L , and is denoted as "lower beam branch" (LBB); a stable branch determined by the local GAM frequency, denoted as the "GAM branch" (GB), and is little affected by the EP distribution function; and a marginally stable branch determined by the birth energy of the distribution function (E b ), and is denoted as "upper beam branch" (UBB). It is noteworthy that, the unstable LBB can have a frequency much higher than the local GAM frequency, while γ/ω G ∝ √ N b can be obtained by balancing the thermal plasma contribution to the dominant destabilizing term due to low energy end cutoff (i.e., the term proportional to 1/(1 − ω 2 L /ω 2 ) in equation (12)). The stability of the three branches described by the dispersion relation (12) on ω L is briefly summarized in Appendix A for the convenience of readers. One important feature is that, the linear unstable LBB frequency is non-perturbatively determined by ω L , and thus, will self-consistently chirp up or down, due to EGAM induced pitch angle scattering, as we show in Appendix B.
Thus, the primary mode in LHD experiment [23] corresponds to the linearly unstable LBB, while the secondary mode with the frequency being local GAM frequency, corresponds to the linearly stable GB. We will show in the next section that, as the unstable LBB chirping up to twice local GAM frequency, it can decay into two linearly stable GBs, similar to the well-known "two plasmon decay" process describing an electromagnetic wave decay into two Langmuir waves [42,43].

Two plasmon decay of EGAM
In this section, the high frequency LBB decay into two linearly stable GBs, will be investigated using nonlinear gyrokinetic theory. Denoting the "primary mode" and two "secondary modes" with subscripts "0", "1" and "2", respectively, the nonlinear particle responses to GAM/EGAM can be derived from nonlinear gyrokinetic equation To properly assess the nonlinear coupling, the contribution of both thermal plasmas and EPs are considered simultaneously, by deriving their nonlinear response to GAM/EGAM, and evaluating their contribution to the surface averaged quasineutrality condition, equation (4).

Negligible thermal plasma contribution to nonlinear coupling
For electrons with ω tr,e ≫ ω G , the nonlinear electron response is dominated by surface averaged contribution. Noting δH L e = (e/T e )F 0,e δφ k , one has i.e., δH N L k,e = 0 up to the O(ω G /ω tr,e ) ≪ 1 order. Thus, electron contribution to the nonlinear coupling can be neglected.
Nonlinear ion response can be derived, noting the ω G ≫ ω tr,i , ω d ordering and that δH L i ≃ −(e/T i )J 0 δφ G (1 + ω d /ω), and one has Note that, for GAM/EAM with n = 0, the geodesic curvature induced drift ω d ∝ sin θ, which gives δH N L k,i a ∝ cos θ dependence to the leading order. Thus, finite contribution to the surface averaged quasi-neutrality can only enter through toroidal coupling (B ∝ 1 − ǫ cos θ), as was discussed in Ref. [50]. Thus, the thermal ion induced nonlinearity via surfaced averaged quasi-neutrality condition, can be estimated to be of order ∼ cǫn 0 k r δφ 0 δφ 1 ω d,i /(ω G rB 0 ), which is comparable to parallel nonlinearity and is, thus, negligibly small. Consequently, we expect the EP induced nonlinear coupling will be dominant, as shown in Refs. [47,48].

Finite EP contribution to nonlinear coupling via EP FoW effects
Nonlinear EP response to the sideband Ω k can be derived from nonlinear gyrokinetic equation, by drift orbit center coordinate transformation. For the Ω 1 generation due to the coupling of Ω 0 and Ω 2 * coupling, taking δH N L 1 = e iΛ1 δH N L d1 , the corresponding equation for nonlinear EP response to Ω 1 can be written as Here, one expects that the contribution from δφ 2 * δH L 0 being dominant with respect to δφ 0 δH L 2 * due to the crucial contribution of resonant EPs (note that Ω 0 is linearly unstable due to resonant EP drive). However, both terms are consistently kept for now, with both resonant and non-resonant EP contribution to the nonlinear coupling accounted for on the same footing. As we will show a posteriori, the nonlinear coupling is dominated by EP FoW effects, while the sub-dominant EP FLR effects can be neglected systematically. Substituting the linear EP response into equation (16), we then have, The nonlinear EP response to Ω 1 can then be obtained, by the pull-back transformation, and one has In the above expression, the first term in the bracket comes from δφ 2 * δH L 0 , as evident from the denominator ω 0 − lω tr , while the other term comes from δφ 0 δH L 2 * . For simplicity, in the following derivation, only the first term due to δφ 2 * δH L 0 will be kept, which can be dominant due to resonant EP contribution. The other term, can also contribute and quantitatively impact the nonlinear process, but we expect its contribution is relatively small. The physics meaning of various terms in equation (17) is clear, in that ω 0 − lω tr in the denominator gives wave-particle power exchanges with the pump Ω 0 , (l+p) comes from ∂ θ δH L 0 in the perpendicular nonlinearity, and the Bessel functions are from EP FoW effects, determining the strength of EPs interaction with the mode.
For finite Ω 1 generation due to δH N L 1 , assuming |Λ k | ≪ 1 for EGAMs with typically global mode structure, one requires the following "selection rules" for nonvanishing EP contribution: 1. l + p + η + h = 0 for non-vanishing component surface average, 2. l + p = 0 for ∂ θ δH L 0 = 0, 3. l = 0 for finite EP linear drive to the primary EGAM, and 4. |l| + |p| + |η| + |h| as small as possible for maximized contribution under theΛ k ≪ 1 assumption. One then has, the dominant contribution comes from l = 1, p = 0, η = −1, h = 0 or l = 1, p = 0, η = 0, h = −1, which gives: The ratio of thermal ion and EP contribution to the nonlinear coupling, obtained from equations (15) and (18), respectively, can be estimated to be confirming the conjecture following equation (15). Here, N i and N h are, respectively, the thermal ion and EP contribution to nonlinear coupling, respectively. Linear EGAM orderings, including β h β i , |γ/ω G | ∼ √ N b and ω tr,h ∼ ω G [26] are used in estimating the ordering of Substituting equation (18) into quasi-neutrality condition, one obtains the equation for Ω 1 generation due to Ω 0 and Ω 2 * coupling Here, E 1 ≡ E EGAM (ω 1 , k r,1 ) is the linear EGAM dispersion relation derived in Ref. [30], which is linearly stable for the GB Ω 1 , with the frequency determined by GAM frequency. Furthermore, A ≡ cΩ 2 i /(rB 0 n 0 ), and The equation for Ω 2 can be similarly derived as with E 2 ≡ E EGAM (ω 2 * , k r,2 * ) being the dispersion relation of Ω 2 * , and The EGAM two plasmon decay dispersion relation can then be derived from equations (19) and (20) as For the two stable GBs, we have, and similarly, with ∆ ≡ ω 1 − ω 0 /2 being the mismatch of half primary mode to GBs, and γ being the secondary mode growth rate due to pump Ω 0 drive. Substituting equations (22) and (23) into (21), we have with the right hand side being the nonlinear EP drive. The secondary modes can be driven unstable as the nonlinear drive overcomes the threshold due to mismatch between half primary mode frequency and GB. Thus, the secondary modes can be excited when the primary amplitude is large enough, or when the frequency mismatch is sufficiently small, i.e., the secondary mode excitation condition is optimized as the primary mode frequency up sweeping to twice local GAM frequency, as observed in the experiment and numerical simulation [37,39].
If we consider only resonant EP contribution to the nonlinear coupling, taking ω tr,h ≃ ω 0 , and assume small EP drift orbit by taking J 0 (Λ) ≃ 1, J 1 (Λ) ≃Λ/2, equation (24) can be reduced to It is clear from equation (25) that, the EPs nonlinearly "drive" the secondary modes are the same particles that resonantly drive the primary mode unstable, though the secondary mode frequency is very different from that of the primary mode. Thus, theoretical results from equation (25) illuminate experimental observations as well as the findings from numerical simulations. For more quantitative comparison, the threshold of primary EGAM amplitude for secondary mode generation can be estimated by with γ G being the GAM collisionless damping rate, and Max(∆, γ G ) giving the maximum value of ∆ and γ G .

Conclusions and Discussions
In this work, an analytical theory is proposed to interpret the secondary mode generation during primary EGAM frequency chirping observed in LHD experiments [37], which is re-produced in MEGA simulation [39]. The interpretation is based on a previous theory on linear EGAM stability in the same LHD low collisionality plasmas, which shows that for the not fully slowed down EP distribution function, the unstable branch (LBB) frequency can be significantly higher than the local GAM frequency; while there is a linearly stable branch (GB) with the frequency determined by the local GAM frequency. The "primary" and "secondary" modes in the experimental observations [23] correspond to the linearly unstable LBB and linearly stable GB, respectively. It is shown that, the LBB can decay into to two linearly stable GBs as its frequency is up-chirping to twice GB frequency, in a process similar to the well known two plasmon decay process [42,43]. The contribution of both thermal plasmas and EPs to the nonlinear process are derived and evaluated. It is found that the thermal plasma contribution to the coupling is negligible compared to that of EP FoWs, and this explains that the nonlinear coupling still occurs when fluid nonlinearity is turned off in the hybrid MEGA simulation. The resonant EPs play crucial roles in the nonlinear coupling, consistent with the observation that the EPs, which "drive" the secondary mode, are the same as those linearly driving the primary mode unstable, though the secondary mode has a frequency much lower than that of the primary. It is noteworthy that the GB frequency is determined by local GAM frequency, and is only slightly modified when GB and LBB strongly couple. Thus, as the primary frequency keeps chirping up, the secondary mode frequency is almost unchanged, as shown in both LHD experimental observation and MEGA simulation. Thus, the present theory, illuminates all the crucial evidence from the experiment and simulation, suggesting this as the mechanism controlling the underlying the physics.
The present theory is local, facilitated by the existence of the high frequency LBB in the low collisionality condition of this specific LHD experiment. However, two plasmon decay process can also occur in typical discharges of usual magnetically confined toroidal plasmas, where the unstable EGAM frequency is typically lower than local GAM frequency. This global two plasmon decay process can occur due to the GAM frequency dependence on local plasma parameters, such that EGAM frequency can be significantly higher than GAM continuum frequency of a distant region, if the temperature gradient is sharp enough. Thus, the two plasmon decay process can occur as the EGAM tunnels through the potential barrier, and strongly couple to GAM where the GAM frequency is half of the EGAM. The condition for the above process to happen is, though, quite difficult to satisfy, since 1. the thermal plasma temperature gradient needs to be sharp, such that the potential barrier is narrow enough to have finite EGAM tunneling through, and 2. the characteristic scale length of thermal plasma temperature nonuniformity needs to be larger than EP density scale length to have EGAM localized by the potential well. The in-depth discussion of the global coupling process is beyond the scope of the present paper focusing on the specific condition of LHD experiments, and will be presented in a separated work.
controlled by three dominant parameters, i.e., ω b and ω L determined by NBI birth energy E b and low energy end cutoff E L , and local GAM frequency ω G . Here, equation (12) was solved for the EGAM stability v.s. ω L , given ω b = 3ω G and λ 0 B 0 < 2/5, such that only the simple pole like singularity is destabilizing. Besides, a small but finite GAM damping rate is assumed γ G = −0.05ω G . Note that ω L = ω b corresponds to all the EPs having the same energy; i.e., to beam ions having not slowed down at all; while ω L ≪ ω b (more precisely, E b > E L ≫ E crit ) corresponds to NBI being fully slowed down. The dependence of the real frequencies of the three branches on ω L is given in Fig. A1, while their growth rates are given in Fig. A2, with the frequencies/growth rates normalized with ω G . It is shown that, UBB frequency is determined by ω b while it remains almost independent of ω G or ω L , and it is marginally stable. The LBB frequency is determined by ω L , and it can be unstable even when its frequency is significantly higher than local GAM frequency. The GB frequency is determined by local GAM frequency, and is almost independent of ω L .
The important information here is that, 1. LBB frequency is determined by ω L ≡ 2E L (1 − λ 0 B 0 )/(qR 0 ) and the unstable LBB frequency can be significantly higher than local GAM frequency if strongly driven; and 2. GB frequency is determined by local GAM frequency, and its dependence on ω L is very weak. These two points are crucial for the interpretation of the nonlinear process in Sec. 4.
Appendix B. Self-consistent EGAM frequency chirping due to pitch angle scattering The nonlinear evolution of EGAM including frequency chirping, can be derived by self-consistently including the slowly evolution of the "equilibrium" EP distribution function F 0,h due to the nonlinear interactions with EGAM [9,10], as addressed in Ref. [18]. The corresponding F 0,h evolution obeys the following Dyson equation [9,10,51,52] ωF 0,h = − e 2ω d 16 Here,ω denotes the slow nonlinear evolution of F 0,h from its initial value F 0,h (0),F 0,h is the Laplace transform of F 0,h , and ω 0,R is the real frequency of EGAM. Equation (B.1) describes the self-consistent evolution of F 0,h , due to emission and re-absorption of a single coherent EGAM. In deriving equation (B.1), only evolution in E (or v or λ) is considered, since both P φ and µ are well conserved for GAM/EGAM with n = 0 and frequency significantly lower than Ω i . Equation (B.1) is derived assuming well circulating EPs andΛ k ≪ 1, while a more systematic treatment can be done following Ref. [10].
EGAM may scatter EPs to smaller pitch angle, and thus, larger ω tr,h regime, which will lead to self-consistent EGAM frequency up-chirping as its frequency is nonperturbatively determined by ω L , as we shown in Appendix A. However, the self-consistent EGAM evolution will be only qualitatively but not quantitatively investigated in the present work, since it will be addressed in an independent work.