Evidence of “Two Plasmon” Decay of Energetic Particle Induced Geodesic Acoustic Mode

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\simeq 0,\pm 1\cdots$, 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\simeq 0,\pm 1,\cdots$, and an up-down anti-symmetric $m\simeq 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 [15, 14, 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\sim\sqrt{L_{n}\rho_{h}}$, with $L_{n}$ being the EP density profile scale length, and $\rho_{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 $\sim$ 7keV) and low plasma density (line averaged density $\sim 0.1\times 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_{\parallel}$ 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 A, while the frequency up-chirping of the primary EGAM, being an important condition for the decay process, is briefly outlined in B.

In the above described nonlinear decay process, a primary EGAM ($\Omega_{0}\equiv\Omega_{0}(\omega_{0},k_{r,0})$) decays into two secondary modes with almost identical frequency, $\Omega_{1}\equiv\Omega_{1}(\omega_{1},k_{r,1})$ and $\Omega_{2}\equiv\Omega_{2}(\omega_{2},k_{r,2})$. All the three modes, are EGAM/GAMs with $n=0/m\simeq 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:

$\displaystyle\delta n_{e}=\delta n_{i}+\delta n_{h},$

(1)

with the perturbed particle density derived from the perturbed distribution function as

$\displaystyle\delta n_{s}\equiv\left\langle(q/m)\partial_{E}F_{0,s}\delta\phi+J_{0}(k_{\perp}\rho_{L,s})\delta H_{s}\right\rangle.$

(2)

Here, subscripts $s=e,i,h$ denote electrons, ions and energetic (hot) particles, respectively, $F_{0,s}$ is the equilibrium distribution function, $\partial_{E}F_{0,s}$ is the short notation for $\partial F_{0,s}/\partial E$, with $E\equiv v^{2}/2$, $J_{0}(k_{\perp}\rho_{L,s})$ is the Bessel function of zero index accounting for finite Larmor radius (FLR) effects, $k_{\perp}$ is the perpendicular wavenumber, $\rho_{L}\equiv v_{\perp}/\Omega_{s}$ is the Larmor radius with $\Omega_{s}$ being gyro-frequency, and $\langle\cdots\rangle$ denotes velocity space integration. $\delta H$ is the non-adiabatic particle response to GAM/EGAM, and can be derived from nonlinear gyrokinetic equation [44]:

	$\displaystyle\left(\partial_{t}+\omega_{tr}\partial_{\theta}+i\omega_{d}\right)\delta H_{k}$	$\displaystyle=$	$\displaystyle-i\omega\frac{q}{m}\partial_{E}F_{0}J_{k}\delta\phi_{k}$		(3)
			$\displaystyle-\frac{c}{B_{0}}\sum_{\mathbf{k}=\mathbf{k^{\prime}}+\mathbf{k^{\prime\prime}}}\hat{\mathbf{b}}\cdot\mathbf{k}^{\prime\prime}\times\mathbf{k}^{\prime}J_{k^{\prime}}\delta\phi_{k^{\prime}}\delta H_{k^{\prime\prime}}.$

Here, $\omega_{tr}\equiv v_{\parallel}/(qR_{0})$ is the transit frequency, $\omega_{d}=k_{r}v_{d}$ is the magnetic drift frequency associated with geodesic curvature, and $v_{d}=(v^{2}_{\perp}+2v^{2}_{\parallel})\sin\theta/(2\Omega R_{0})\equiv\hat{v}_{d}\sin\theta$ is the magnetic drift velocity, $J_{k}\equiv J_{0}(k_{\perp}\rho_{L})$ for simplicity of notation, and $\sum_{\mathbf{k}=\mathbf{k^{\prime}}+\mathbf{k^{\prime\prime}}}$ 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) reduces to

$\displaystyle-\frac{e}{m_{i}}n_{0}k^{2}_{r}\frac{1}{\Omega^{2}_{i}}\left(1-\frac{\omega^{2}_{G}}{\omega^{2}}\right)\overline{\delta\phi}_{k}+\overline{\delta n}^{NL}_{i}+\overline{\delta n}_{h}=0,$

(4)

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 $\overline{(\cdots)}\equiv\int^{2\pi}_{0}(\cdots)d\theta/(2\pi)$ is the surface averaging. Furthermore, $\omega_{G}\equiv\sqrt{7/4+\tau}v_{it}/R$ is the leading order GAM frequency with $\tau=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 $\omega_{G}$ with more accurate expression [20]. The linear EGAM dispersion relation can be derived by keeping $\delta 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, 25, 26, 30, 29]. Nonlinear modulation of GAM/EGAM by thermal plasma nonlinearity, can be accounted for by $\delta n^{NL}_{i}$, including excitation by DW/DAWs [14, 45, 16]. 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 [48, 47]. The thermal plasma contribution to the nonlinear coupling [46], meanwhile, will be shown to be negligible compared to the dominant role of EPs.

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 $\delta H_{k}=\delta H^{L}_{k}+\delta H^{NL}_{k}$, the linear EP response to EGAM can be derived from the linear gyrokinetic equation,

$\displaystyle\left(\partial_{t}+\omega_{tr}\partial_{\theta}+i\omega_{d}\right)\delta H^{L}_{k}=-i\omega\frac{e}{m}\partial_{E}F_{0}J_{k}\delta\phi_{k}.$

(5)

Equation (5) can be solved by transforming into the EP drift orbit center coordinate by taking $\delta H^{L}_{k}\equiv e^{i\Lambda_{k}}\delta H^{L}_{dk}$, with $\Lambda_{k}$ satisfying $\omega_{tr}\partial_{\theta}\Lambda_{k}+\omega_{d,k}=0$. Here, for simplicity of discussion and focus on proof of principle demonstration, well circulating EPs are assumed, thus, variation of $v_{\parallel}$ along the magnetic field is neglected, and one has $\Lambda_{k}=\hat{\Lambda}_{k}\cos\theta$, with $\hat{\Lambda}_{k}=k_{r}\hat{\rho}_{d}$ being radial wave-number normalized to drift orbit width, and $\hat{\rho}_{d}=\hat{v}_{dr}/\omega_{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, $\tau\equiv T_{e}/T_{i}\ll 1$ is assumed for simplicity, such that one has $\delta\phi_{G}\simeq\overline{\delta\phi}_{G}$ while $\omega_{tr,e}\gg\omega_{G}$ is still satisfied ¹¹1Interested readers may refer to Ref. [20] for the contribution of finite $\tau$ to linear GAM dispersion relation via the $m\neq 0$ components of the scalar potential.. We then have,

$\displaystyle\left(\partial_{t}+\omega_{tr}\partial_{\theta}\right)\delta H^{L}_{dk}=-\frac{e}{m}e^{-i\hat{\Lambda}_{k}\cos\theta}J_{k}\partial_{t}\overline{\delta\phi}_{k}\partial_{E}F_{0,h},$

(6)

from which $\delta H^{L}_{dk}$ can be derived as

$\displaystyle\delta H^{L}_{dk}=-\frac{e}{m}\partial_{E}F_{0,h}J_{k}\overline{\delta\phi}_{k}\sum^{\infty}_{l=-\infty}\frac{\omega}{\omega-l\omega_{tr}}(-i)^{l}J_{l}(\hat{\Lambda}_{k})e^{il\theta}.$

(7)

Here, $l$ is integer, and $\sum_{l}\equiv\sum_{l=-\infty}^{\infty}$ will be used later for simplicity of notation. In deriving equation (7), the Jacobi-Anger expansion $e^{i\hat{\Lambda}_{k}\cos\theta}=\sum_{l}i^{l}J_{l}(\hat{\Lambda}_{k})e^{il\theta}$ is used. Pulling back to the EP guiding center coordinate, we then have, the linear well-circulating EP response to EGAM

$\displaystyle\delta H^{L}_{k}=-\frac{e}{m}\partial_{E}F_{0,h}J_{k}\overline{\delta\phi}_{k}\sum_{l,p}\frac{\omega}{\omega-l\omega_{tr}}i^{-l+p}J_{l}(\hat{\Lambda}_{k})J_{p}(\hat{\Lambda}_{k})e^{i(l+p)\theta}.$

(8)

Substituting equation (8) into the surface averaged quasi-neutrality condition, equation (4), one then obtains,

$\displaystyle\frac{e}{m_{i}}n_{0}k^{2}_{r}\frac{1}{\Omega^{2}_{i}}\mathscr{E}_{EGAM}\overline{\delta\phi}=0,$

(9)

with the EGAM dispersion relation given by

	$\displaystyle\mathscr{E}_{EGAM}\equiv$		$\displaystyle-1+\frac{\omega^{2}_{G}}{\omega^{2}}$		(10)
			$\displaystyle+\frac{\sqrt{2}\pi B_{0}e^{2}}{n_{0}}\int\frac{(2-\lambda B_{0})^{2}}{(1-\lambda B_{0})^{1/2}}\frac{E^{5/2}dEd\lambda\partial_{E}F_{0,h}}{2E(1-\lambda B_{0})-\omega^{2}q^{2}R^{2}_{0}}.$

Here, only the $l=\pm 1$ transit harmonic are kept in equation (10), in consistency with the typical $\hat{\Lambda}_{k}\ll 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

$\displaystyle F_{0,h}=\frac{c_{0}H(E_{b}-E)H(E-E_{L})}{E^{3/2}+E^{3/2}_{crit}}\delta(\lambda-\lambda_{0}),$

(11)

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, $\lambda\equiv\mu B_{0}/E$ is the pitch angle, $\mu\equiv mv^{2}_{\perp}/B$ is the magnetic moment, $H(\lambda-\lambda_{0})$ is the Heaviside step function, $E_{b}$ is the EP birth energy, $E_{L}\simeq E_{b}\exp(-2\gamma_{c}t)$ 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 $\gamma_{c}$ [49]. The normalization to EP density $n_{b}$ gives $c_{0}=n_{b}\sqrt{1-\lambda_{0}B_{0}}/(2\sqrt{2}\pi B_{0}\ln(E_{b}/E_{L}))$. Note that, $c_{0}$ is proportional to $\Gamma/\gamma_{c}$ with $\Gamma$ being the NBI particle flux and $\gamma_{c}$ is the slowing-down rate on thermal electrons. The resulting dispersion relation is given as [30]

	$\displaystyle-1+\frac{\omega^{2}_{G}}{\omega^{2}}+N_{b}\left[C\left(\ln\left(1-\frac{\omega^{2}_{b}}{\omega^{2}}\right)-\ln\left(1-\frac{\omega^{2}_{L}}{\omega^{2}}\right)\right)\right.$
	$\displaystyle\left.+D\left(\frac{1}{1-\omega^{2}_{b}/\omega^{2}}-\frac{1}{1-\omega^{2}_{L}/\omega^{2}}\right)\right]=0.$		(12)

Here, $N_{b}\equiv n_{b}q^{2}\sqrt{1-\lambda_{0}B_{0}}/(4\ln(E_{b}/E_{L})n_{0})$ is the ratio of EP to bulk plasma density, $C=(2-\lambda_{0}B_{0})(5\lambda_{0}B_{0}-2)/(2(1-\lambda_{0}B_{0})^{5/2})$, $D=\lambda_{0}B_{0}(2-\lambda_{0}B_{0})^{2}/(1-\lambda_{0}B_{0})^{5/2}$, and $\omega_{b}\equiv\sqrt{2E_{b}(1-\lambda_{0}B_{0})}/(qR_{0})$ and $\omega_{L}\equiv\sqrt{2E_{L}(1-\lambda_{0}B_{0})}/(qR_{0})$ 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 $\lambda_{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 $\omega_{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 $\gamma/\omega_{G}\propto\sqrt{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-\omega^{2}_{L}/\omega^{2})$ in equation (12)). The stability of the three branches described by the dispersion relation (12) on $\omega_{L}$ is briefly summarized in A for the convenience of readers. One important feature is that, the linear unstable LBB frequency is non-perturbatively determined by $\omega_{L}$, and thus, will self-consistently chirp up or down, due to EGAM induced pitch angle scattering, as we show in 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].

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

$\displaystyle\left(\partial_{t}+\omega_{tr}\partial_{\theta}+i\omega_{d}\right)\delta H^{NL}_{k}=-\frac{c}{B_{0}}\sum_{\mathbf{k}=\mathbf{k}^{\prime\prime}+\mathbf{k}^{\prime}}\hat{\mathbf{b}}\cdot\mathbf{k}^{\prime\prime}\times\mathbf{k^{\prime}}J_{k^{\prime}}\overline{\delta\phi}_{k^{\prime}}\delta H_{k^{\prime\prime}}.$

(13)

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 quasi-neutrality condition, equation (4).

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.

This work is supported by the National Key R&D Program of China under Grant No. 2017YFE0301900, the National Science Foundation of China under grant No. 11875233. This work was also carried out within the framework of the EUROfusion Consortium and received funding from the EURATOM research and training programme 2014 - 2018 and 2019 - 2020 under Grant Agreement No. 633053 (Project Nos. WP19-ER/ENEA-05 and WP17-ER/MPG-01). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The data that support the findings of this study are available from the corresponding author upon reasonable request.