Global simulations of kinetic-magnetohydrodynamic processes with energetic electrons in tokamak plasmas

The energetic electrons (EEs) produced from auxiliary heatings have been found to destabilize various Alfven eigenmodes (AEs) in recent experiments. To investigate EE relevant kinetic-magnetohydrodynamic (MHD) processes, a global fluid-kinetic hybrid model is formulated and verified in this work, which consists of Landau-fluid bulk plasmas and drift-kinetic EEs that incorporates the dominant damping and drive mechanisms, respectively. The numerical capability of Landau-fluid bulk plasmas is obtained based on a well-benchmarked eigenvalue code multiscale analysis for plasma stability (MAS) using general geometry (Bao et al 2023 Nucl. Fusion 63 076021), and the comprehensive EE responses to the low frequency ( ω≪Ωci ) MHD fluctuations are analytically derived and implemented in MAS, which not only cover contributions from trapped and passing particles, but also take into account the effects of adiabatic fluid convection and non-adiabatic kinetic compression. The linear properties of EE-driven beta-induced AEs (e-BAEs) are systematically studied using both MAS model and gyrokinetic toroidal code (GTC) particle-in-cell simulations. Building on the good agreements on the mode structure and dispersion relation, several key issues of e-BAE physics are analyzed and discussed, including the parametric dependences of e-BAE stability on EE mass, temperature and density with corresponding phase space dynamics, EE non-perturbative effects on the symmetry breaking of mode structure, and the EE density and temperature thresholds for e-BAE excitation that overcome bulk plasma damping. With these efforts, the upgraded MAS model is superior than initial-value simulations restricted by stringent electron Courant condition for fast linear analyses of most EE-AE problems with ω≪Ωci , of which wave–particle resonance can be used to analyze the analogous effect of alpha particle driven AEs in future fusion reactor.


Introduction
Various Alfven eigenmodes can be destabilized by energetic particles (EPs) through abundant channels of waveparticle resonances in tokamak plasmas, and the relevant kinetic-magnetohydrodynamic (MHD) processes are of great interest from aspects of experiment, theory and simulation [1,2].In particular, the AEs driven by energetic ions (EIs) are widely observed in experiments, and most linear and nonlinear properties have been studied and understood over the past four decades [3,4], while the AEs driven by energetic electrons (EEs) are much less explored before until recently, i.e., not only the single mode dynamics of EE-driven betainduced Alfven eigenmode (e-BAE) and EE-driven toroidal Alfven eigenmode (e-TAE) are diagnosed in HL-2A [5,6] and in EAST [7,8,9,10], but also the experimental evidences of the nonlinear mode-mode interaction between e-TAE and geodesic acoustic mode (GAM) have been demonstrated [11].Understanding the excitation mechanism and required plasma condition of EE-driven AEs are helpful for clarifying the relevant physics in these present-day experiments, moreover, the results can be analogous to the alpha particle physics in future fusion power plant characterized by the similarly small dimensionless orbits.
Both theoretical and numerical efforts have been made to explain the EE-driven AEs in experiments.It has been shown that the bounce-averaged dynamics of EEs are responsible for the resonant interaction with MHD modes characterized by much lower frequency than EE bounce and transit frequencies (ω ω b,h < ω t,h ) from the perspective of first-principle gyrokinetic framework [12].Recently, the bounce-kinetic EE model is also applied to the theoretical studies of linear excitation of e-BAE [13,14] and the nonlinear excitation of zonal flow by pump e-BAE [15], which successfully reveal that the dominant destabilizing mechanism of e-BAE is attributed to the EE precessional drift resonance, as well as the important role of resonant EEs on the nonlinear mode-mode coupling.On the other hand, both the gyrokinetic particle-in-cell (PIC) simulation [16,17] and kinetic-MHD hybrid simulation [18,19] have been performed to study the e-BAE and e-TAE with drift-kinetic description of EE dynamics, and the deeply-trapped EE response is found to have a maximal resonance island with e-BAE [16] and e-TAE [18] propagating along the electron diamagnetic direction, while the passing EEs are responsible for the high frequency EAE [18].However, regarding to the shot-to-shot analysis in experiments, the ballooning theory relies on the scale separation and is difficult to incorporate the complex geometry, and the initial value simulation using drift-kinetic EEs suffers the stringent and unnecessary numerical constraints due to the electron Courant condition with realistic electron-ion mass ratio [20].For comparison, the eigenvalue approach using magnetic coordinates is more efficient for the plasma stability analysis with realistic geometry, such as LIGKA [21], NOVA-K [22] and CASTOR-K [23], which mostly focus on EI physics rather than EE physics.Motivated by previous theories and simulations that have provided the fundamental physics insights on the linear excitation of AEs by EE preccessional drift resonance in experiments, we construct kinetic-MHD model associated with EE physics using eigenvalue approach in this work, which combines the necessary EE drive, bulk plasma damping and realistic geometry together for analyzing and optimizing experiments with fast parameter scans.
MAS (Multiscale Analysis for plasma Stabilities) eigenvalue code has been developed for plasma stability analysis in experimental geometry based on a five-field Landau-fluid model description of bulk plasma, which incorporates kinetic effects into MHD modes including ion finite Larmor radius (FLR), ion and electron diamagnetic drifts and Landau dampings on the same footing with fluid responses [24].In this paper, we further extend MAS formulation to include EE physics described by bounce-averaged drift-kinetic equation, and a modified deeply-trapped model is proposed to solve the EE perturbed distribution and calculate corresponding kinetic particle compression (KPC) term, which is then coupled to the existing numerical capability of bulk plasma in a non-perturbative manner.After benchmarking with first-principle gyrokinetic simulations, the upgraded MAS code is useful for explaining EE-driven AEs in relevant experiments as well as providing quasilinear spectra for integrated modeling such as machine learning [25].The remainder of this paper is organized as follows.The formulation of hybrid model that couples drift-kinetic EE and Landau-fluid bulk plasma is introduced in section 2. The EE responses in various polarization and frequency limits of electromagnetic fields are shown in section 3.In section 4, the validity regime of deeply-trapped approximation is delineated, and the verification of the implemented EE module in e-BAE simulations is presented.The conclusions are given in section 5.

Drift-kinetic equation for energetic electron
Considering the smallness of realistic electron mass, we apply the drift-kinetic model for EE species that induces the necessary Landau resonance while ignores the FLR effect.The linearized drift-kinetic equation in five dimensional phase space uses guiding center position R, magnetic moment µ = m e v 2 ⊥ /2B 0 and parallel velocity v || as independent variables, which reads where δf h and f h0 are the EE perturbed and equilibrium distributions, L 0 and δL L are the equilibrium and linear perturbed propagators given by q e and m e are electron charge and mass, δφ and δA || are the electrostatic potential and parallel vector potential, b 0 = B 0 /B 0 is the unit vector of the equilibrium magnetic field, 0 /cm e is the electron cyclotron frequency, δB = ∇ × δA || b 0 represents the perturbed magnetic field, and v d and v E are the magnetic drift and E × B drift, which read For EE-driven Alfvenic instabilities, the linear unstable spectra can be influenced by the EE velocity distribution relying on the specific auxiliary heating methods such as ECRH and LHCD [12].Though the EE velocity distribution is not unique in experiments, we define effective EE temperatures for analysis convenience, namely, || f h0 dv/n h0 and T ⊥h0 = µB 0 f h0 dv/n h0 , where n h0 = f h0 dv is the EE density.In current work that illustrates the model scheme, we utilize the isotropic Maxwellian to approximate the EE equilibrium distribution, i.e., with T h0 = T ||h0 = T ⊥h0 , and other EE velocity distributions such as slowing-down will be investigated in future work.We further separate the perturbed distribution δf h in Eq. ( 1) into the adiabatic part and non-adiabatic part The adiabatic distribution is determined by the terms related to fast parallel dynamics in Eq. ( 1), which are in the leading order of ω/(k || v || ) as where For analysis convenience, we newly define δψ in terms of δA || through relation Substituting Eq. ( 4) into Eq.( 3) and using the ansatz ∂ t = −iω and b 0 • ∇ = ik || , δf A can be readily solved from Eq. ( 3) where Note that the form of δf A in Eq. ( 5) is also adopted in gyrokinetic theory [26,27] and simulation [28], which contains both the adiabatic response to the parallel electric field and convective response to the perturbed magnetic field.
From Eqs. (1), ( 2) and ( 5), the governing equation for the non-adiabatic distribution δK can be written as where where we have δK p ∼ 0 since v || b 0 • ∇ → ∞, and δK p contribution to perturbed density and pressure can be ignored.However, v || δK p is finite and gives rise to a fraction of parallel current density carried by passing EEs.
Second, EE particles interact and exchange energy with MHD modes primarily through the precessional drift resonance, while the bounce and transit frequencies are too high to resonate with typical AEs, as demonstrated and confirmed by recent theories and simulations [13,15,16,18].According to these characteristics, we have ω ∼  [29], and then perform bounce-average indicates that δK t δK t b and finite orbit width (FOW) can be dropped for EE particles as a higher order effect, then the solution of trapped EE non-adiabatic response is where term {I} corresponds to the non-ideal MHD effect ∆φ = δφ − δψ, and term {II} drives the badcurvature instabilities with a minimum threshold ω * p,h > ω.For most instabilities in tokamak characterized with 'flute-like' mode structures, δφ and δψ peak around the rational surfaces where |nq − m| 1, then we can further simplify Eq. ( 8) using δφ ≈ δφ, δψ ≈ δψ, and ω d δψ ≈ ω d δψ.The more general validity regime of these simplifications is exp [−i (m − nq) θ] ≈ 1 [15], which only requires |θ| 1/|nq − m|.Thus, for instabilities with moderate |nq − m| ≤ 1 that deviate from rational surface, one can reduce trapped EE fraction by decreasing the θ domain of banana orbit mirror throat to more deeply-trapped regime, which still hold δφ ≈ δφ, δψ ≈ δψ, and ω d δψ ≈ ω d δψ for model accuracy.(Note that these simplifications are applied for deriving the EE moments in section 2.2).

Energetic electron continuity equation
To couple drift-kinetic EEs with the existing Landau-fluid model of bulk plasmas, one needs to derive the EE continuity equation by integrating Eq. ( 1) in velocity space using where and Taking the moment of Eq. ( 5), the adiabatic components in Eqs. ( 10) and ( 11) are obtained as and Eqs. ( 12) -( 15) are the adiabatic EE moments which will be used in the hybrid model.Then we shall derive the non-adiabatic EE moments.Since the passing EEs move much faster than Alfven speed of which response is almost adiabatic to shear Alfven wave (SAW) and ion sound wave (ISW), δK p in Eq. ( 7) does not contribute to density and pressure perturbations but can lead to a finite current in high-v || regime.On the other hand, the precessional drift of trapped EE can be close to or below the Alfven speed, and δK t in Eq. ( 8) are responsible for the non-adiabatic EE density and pressure.Moreover, it has been shown that the deeply-trapped EEs can effectively interact with most drift-type instabilities through precessional drift resonance [13,16,18], thus we apply the deeply-trapped approximation for computing the precession frequency in Eq. ( 8) where ω d0 denotes the normal curvature component of ω d = −iv d • ∇ at the outer midplane (θ = 0), while the geodesic curvature component does not contribute to the bounce-averaged quantity.The specific form of ω d0 in flux coordinates is given by Eq. ( 69), and the validation of Eq. ( 16) is discussed in section 4. Meanwhile, with the deeply-trapped approximation, the last term in Eq. ( 9) becomes where ω D0 = T h0 meE ω d0 , and δP N A h 1 2 trap δK t Edv represents the non-adiabatic EE pressure.We use pitch angle λ = µB a / (m e E) and energy per unit mass E = v 2 /2 as the two dimensions (λ, E) of velocity space instead of v || , µ , where B a is the on-axis magnetic field strength, then the velocity space integration with trapped particles can be expanded as where σ = sign(v || ), and B a /B max ≤ λ low ≤ B a /B 0 is the lower cutoff of trapped particle pitch angle.For Maxwellian equilibrium distribution f h0 = n h0 me 2πT h0 3/2 exp − meE T h0 , the trapped particle fraction f t at B 0 location can be expressed as By counting all trapped particles at B 0 location with λ low = B a /B max , it is seen that f t = 0 in the high field side with B 0 = B max and f t = 1 − B min /B max in the low field side with B 0 = B min .Compared to former gyrokinetic theory that applies f t ≈ √ 2 (where = r/R 0 ) and ignores poloidal variation [30], Eq. ( 19) has both radial and poloidal dependences and reflects the realistic trapped particle distribution in the poloidal plane.In general, Eq. ( 18) can be simplified when the integrands do not rely on λ Thus, with the assumption in Eq. ( 16), it is straightforward to integrate δK t in Eq. ( 8) in velocity space using Eq. ( 20), and the non-adiabatic density and pressure perturbations are derived explicitly and where ζ = ω/ω D0 , and the response functions are given by Substituting Eqs. ( 12)-( 15), ( 17), ( 21) and ( 22) into Eq.( 9), the non-adiabatic parallel velocity δu N A ||h is readily solved as where ω D = T h0 meE ω d .In the derivation of Eq. ( 23), the resonance terms associated with Z √ ζ in Eqs. ( 21) and ( 22) exactly cancel out with each other through Eq. ( 9), which again proves that the trapped EE particles do not carry parallel current as we assumed before.On the other hand, δu N A ||h can also be computed by integrating Eq. ( 7) in velocity space with passing EE fraction, which is consistent with Eq. ( 23) in the leading order.

Formulation of Landau-fluid bulk plasma and drift-kinetic energetic electron hybrid model
Next, we couple the EE moments described by Eqs. ( 12)-( 15) and ( 21)-( 23) to the Landau-fluid model of bulk plasma in Ref. [24], and then yield the fluid-kinetic hybrid model as follows where δP A h = δP A ||h = δP A ⊥h in Eq. ( 24), Z i n i0 +q e (n e0 + n h0 ) = 0, and the definitions of bulk plasma variables are consistent with Ref. [24].The two-moment and three-moment Landau closures are applied to thermal electrons in Eq. ( 25) and thermal ions in Eq. ( 26) respectively, which show good agreements with drift-kinetic theory on the response functions in the regime of [31].To close above equation set, we also have the equations for δP e , δT e and δT i as and Meanwhile, the thermal electron perturbed density δn e and parallel velocity δu ||e in Eqs. ( 24), (25) and Eq. ( 27) are calculated through the quasi-neutrality condition and parallel Ampere's law as and Eqs. ( 4), ( 12)-( 15), ( 21)-( 23) and ( 24)-( 33) form a closed system for the fluid-kinetic hybrid simulation model as shown in figure 1.The main characteristics are briefly summarized here, first the well-circulating and deeply-trapped approximations on EE moments do not break the continuity equation, second the EE and bulk plasma are coupled through the quasi-neutrality condition, parallel Ampere's law and vorticity equation in a non-perturbative manner, third the EE-KPC term is responsible for the dissipative excitation of AEs while the EE-IC term contributes to the reactive MHD interchange drive in Eq. ( 24).

Comparison of the moment ordering between energetic electron and thermal electron
To estimate the ordering of each EE moment in the hybrid model, one need to compare δn h , δu ||h and δP h with corresponding thermal electron moments δn e , δu ||e and δP e .Note that Eq. ( 32) for δn e is not intuitive to compare with Eqs. ( 14) and ( 21), we solve δn e by using Eqs.( 4) and (25) equivalently Then the thermal electron pressure δP e is derived using Eqs.( 30) and (34) as where The thermal electron continuity equation can be obtained from Eqs. ( 9), ( 24), ( 28), ( 32) and ( 33) Substituting Eqs. ( 34) and ( 35) into Eq.( 36) and considering ∂ t = −iω and ∇ = ik, we have where 12)-( 14) and ( 23) already have similar forms with Eqs. ( 34), ( 35) and ( 37 and Since EE equilibrium pressure is comparable to thermal electron equilibrium pressure, we can define the smallness parameter and where ω * ,h = ω * n,h + ω * T,h .Moreover, the EE precessional drift resonance requires Then the orders of EE moments can be estimated based on Eqs. ( 42)-(45).In the electrostatic limit δψ = 0, we have and It is seen that only δP N A h is probably close to δP e when ζ → ∞, however, the vorticity Eq. ( 24) that contains δP N A h term becomes redundant in the electrostatic limit, thus EEs are not important to the electrostatic polarized modes.In the electromagnetic limit δφ = δψ, we have and Therefore, δP A h , δP N A h and δu ||h are the leading order terms in the electromagnetic limit, and should be kept for EE species in the hybrid model, while the EE density perturbations δn A h and δn N A h can be safely dropped as higher order small terms.

Verification and validation of the fluid-kinetic hybrid model
The fluid-kinetic hybrid model in section 2.3 can be casted into a nonlinear eigenvalue equation in ω as where A and B are the operators of bulk plasma Landau-fluid model and are independent of ω, and C (ω) is the drift-kinetic EE operator that relies on ω nonlinearly.We use Newton's iterative method to solve Eq. ( 60) with the initial guesses of (ω, X) obtained from AX = ωBX, which is efficient when EE term is small and perturbative.For cases that EE term is comparable with bulk plasma terms, we solve AX−ωBX+ C (ω) X = 0 in the middle steps instead of directly solving Eq. ( 60), where is an adjustable parameter that gradually increases from 0 to 1, thus the proper initial guesses of (ω, X) can be obtained for each -value case and guarantee the robustness of Newton's iterative method in the non-perturbative regime.
To verify the EE physics model and numerical scheme, we first validate the EE precession frequency that uses deeply-trapped approximation in both analytic and experimental tokamak geometries, and then carry out the verification simulations of e-BAE to demonstrate the correctness of EE precessional drift resonance.

How good is deeply-trapped approximation for bounce-averaged guidingcenter dynamics?
As claimed by Eq. ( 16) in section 2, we utilize the deeply-trapped approximation for calculating the precession frequency of trapped EE.To delineate its regime of validity, we compare Eq. ( 16) with the exact precession frequency by performing bounce-average along the realistic banana orbit.The Boozer coordinates (ψ, θ, ζ) is applied in MAS to describe the magnetic field, which has the contravariant form B 0 = q (ψ) ∇ψ×∇θ−∇ψ×∇ζ and the covariant form B 0 = I (ψ) ∇θ + g (ψ) ∇ζ.Then the guiding-center equation of motion in equilibrium B-field can be expressed in (ψ, θ, ζ) coordinates [32] and where α denotes the particle species, J = (∇ψ × ∇θ • ∇ζ) −1 is the Jacobian and ρ || = v || /Ω cα .Given energy E = 0.5v 2 and pitch angle λ = µB a / (m α E) of a particle, parallel velocity v || can be written as and the bounce period can be calculated using Eq. ( 62) and (65) [32] Using Eqs. ( 61) -( 64) and (66) , the exact precession frequency is then derived as [32] where n is the toroidal mode number, and the first term on the RHS is the leading order term that depends on both the radial ψ and poloidal θ coordinates, i.e., ω d = ω d (ψ, θ).
On the other hand, the magnetic drift frequency ω d in Boozer coordinate can be expressed as where term {I}, term {II} and term {III} represent the geodesic curvature, normal curvature and equilibrium current contributions respectively.In Eq. ( 68), we note that term {I} doesn't contribute to the precession frequency since the geodesic drift cancels over a bounce period, effect of term {III} is ignorable, while term {II} is equal to Eq. ( 67) in the limit of τ b → 0, θ → 0 and v || → 0 (i.e., deeply-trapped limit) which is used for the deeply-trapped approximation in Eq. ( 16), and only has ψ dependence, i.e., ω d0 = ω d0 (ψ).
In order to elucidate the accuracy of ω d0 , we compare Eqs. ( 67) and (69) in both the analytic geometry [16] and DIII-D general geometry [33].The safety factor q and magnetic shear s = (1/q)(dq/dr) are given in figure 3.For analysis convenience, we use (ψ b , θ b ) coordinates to represent trapped particles, where ψ b ∈ [0, ψ edge ] is the poloidal magnetic flux of banana orbit center and θ b ∈ [0, 2π] is the Boozer poloidal angle at the banana tip location (i.e., the mirror throat).Using the magnetic field strength at (ψ b , θ b ) location, i.e., B 0,b = B 0 (ψ b , θ b ), we can express the pitch angle of trapped particle as λ = B a /B 0,b , which varies in the range of B a /B max < λ < B a /B min .Substituting Eq. (65) into Eq.( 67), the exact precession frequency ω d can be calculated for all trapped particles in terms of (ψ b , θ b ).The ratio between ω d (ψ b , θ b ) and ω d0 (ψ b ) in the analytic and DIII-D general geometries are shown in figures 4 (a) and 5 (a) respectively.First, it is seen that the numerical calculation of Eq. ( 67) agrees with the measurement from the realistic orbits by Eqs. ( 61)-(64).Second, most areas on the low field sides are characterized with ω d /ω d0 ∼ 1, thus ω d0 in Eq. ( 69) can well approximate the precession frequency for most trapped particles not only in the analytic geometry with low magnetic shear but also in the experimental geometry with broader range of magnetic shear, where the significance of magnetic shear term in Eq. ( 67) is decreased by the elongation effect.Moreover, since the trapped particle fraction f t determines the EE-drive intensity through Eqs. ( 21) and (22), one needs to adjust λ low in Eq. ( 19) so that f t only incorporates the contribution of trapped particles that satisfy ω d /ω d0 ∼ 1, while the barely-trapped particles beyond the model capability should be excluded.Besides the constraint arising from the approximation on precession frequency, we note that the bounce-average operations on electromagnetic fields in Eq. ( 8) are removed for integrating the EE moments, and corresponding validity regime of |θ| 1/|nq − m| becomes another constraint for λ low .Figures 4 (b) and 5 (b) show the 2D profile of f t using λ low = B a /B max in the poloidal plane, which take into account all trapped particles and overestimate the deeply-trapped EE-drive intensity.We then apply λ low = 1 in Eq. ( 19) to calculate f t as shown in figures 4 (c) and 5 (c), which corresponds to the trapped particles with entire banana orbits on the low field side, in consistency with the validity regime of ω d /ω d0 ∼ 1 in figures 4 (a) and 5 (a).Thus, the 2D function f t calculated by using λ low = 1 in Eq. ( 19) correctly reflects the deeply-trapped EE-drive intensity for modes peak at k || ∼ 0 and is applied in the following e-BAE simulations.For modes peak at moderate k || ≤ 1/qR 0 such as e-TAE, one needs to choose λ low < 1 in Eq. ( 19) to calculate f t in the more deeply-trapped regime (i.e., reduce θ b ).For comparison, the widely used simple 1D form of f t ≈ √ 2 assumes trapped particles distribute uniformly along the poloidal direction and might amplify the drive of deeply-trapped fraction.In addition, the omitted barely-trapped EE-drive is in the higher order due to its small population.

Linear properties of beta-induced Alfven eigenmode driven by energetic electrons (e-BAE)
To verify the global fluid-kinetic hybrid simulation model and corresponding numerical implementation in MAS code, we carry out simulations of e-BAE in analytic geometry based on a well-established benchmark case [16], where the safety factor profile and concentric-circular geometry are shown in figures 3 (a) and 4 respectively.Protons are used for thermal ions with Z i = e.The thermal electron, thermal ion and EE temperatures are uniform with T i0 = T e0 = 500eV and T h0 = 25T e0 , and the thermal electron density is uniform with n e0 = 1.3 × 10 14 cm −3 .The EE density profile is described by n h0 = 0.05n e0 1 + 0.2 tanh 0.26 − ψ /0.06 − 1.0 , where ψ = ψ/ψ w is the poloidal magnetic flux normalized by the wall value, and the reciprocal of EE density scale length L −1 n,h = (∇n h0 /n h0 ) −1 peaks at q = 2 rational surface with R 0 /L nh = 12.7.The thermal ion density is then determined by the quasi-neutrality condition Z i n i0 + q e (n e0 + n h0 ) = 0.The on-axis magnetic field strength is B 0 = 1.91T , the major radius is R 0 = 0.65m, the minor radius is a = 0.333R 0 , and the safety factor profile in figure 3 (a) is q = 1.797 + 0.8 ψ − 0.2 ψ2 .The global mode structure and dispersion relation of n = 3 e-BAE are analyzed as follows.The 2D poloidal mode structure of electrostatic potential δφ is shown in figure 6 (a1), which is characterized with a 'boomerang' shape.The dominant principal poloidal harmonic of m = 6 peaks at q = 2 rational surface, of which amplitude is much larger than the neighboring sideband harmonics of m = 5 and m = 7 that corresponds to a weakly ballooning structure as shown in figure 6 (a2).Different from δφ, the poloidal mode structure of parallel vector potential δA || exhibits anti-ballooning feature as shown in figure 6 (b1), where the real parts of m = 5 and m = 7 poloidal sidebands are comparable to the resonant m = 6 harmonic but with phase difference as shown in figure 6 (b2).Moreover, the mode structure of ∆φ = δφ − δψ is shown in figures 6 (c1) and (c2), which represents the derivation from the ideal-MHD limit and reflects the polarization of each poloidal harmonics.The m = 6 harmonic is predominantly Alfvenic according to the small ∆φ since the electrostatic component δφ and electromagnetic component δψ nearly cancel out, while ∆φ perturbation is large in m = 5 and m = 7 sidebands which indicates the acoustic component becomes more important.The symmetry breaking of AE mode structure due to non-ideal MHD effects, such as kinetic EPs induced anti-Hermitian part of dielectric tensor [34], have been widely observed in simulations [35,36] and experiments [37,38], which are characterized with twisted mode structures .To understand the role of EEs on 'boomerang' shape mode structure in figure 6 (a) that is no longer up-down symmetric, four simulation cases are performed with different EE physics as described in table 1. Case {I} corresponds to the Landau-fluid simulation of bulk plasma without any EE effects, and the poloidal mode structure is relatively up-down symmetric as shown in figure 7 (a1), where the small distortion is due to the kinetic effects of bulk plasmas that lead to anti-Hermitian contribution.In case {II}, the EE-KPC term is added on top of case {I}, and the 2D mode structure is radially broadened and drastically twisted with obvious tails on both sides of mode rational surface in figure 7 (b1), which indicates that the anti-Hermitian contribution from EE-KPC is much larger than bulk plasmas.For comparison, figure 7 (c1) shows the 2D mode structure of case {III} that includes EE-IC term on top of case {I} instead of EE-KPC term, and it can be seen that EE-IC term has little impact on the symmetry breaking which is different from case {II}, but EE-IC term can broaden the radial width to a certain extent.In case {IV}, the non-perturbative effects of both EE-IC and EE-KPC terms are considered (figure 7 (d1) and figure 6 (a)), it is seen that the degree of symmetry breaking in case {IV} is between case {II} and case {III}, and the radial width is close to both case {II} and case {III}.Thus, there are two EE non-perturbative effects on e-BAE mode structure: (i) the EE-KPC term provides the dominant kinetic effects and induces the large ant-Hermitian part of dispersion relation, which is responsible for the distortion e-BAE mode structure; (ii) the EE-IC term represents the fluid-like convective response and enhances the Hermitian part, which not only compensates the symmetry breaking induced by EE-KPC term but also broadens the radial width through MHD interchange effect.
On the other hand, the radial profile of phase angle θ r can also reflect the mode structure symmetry, which has received interest from recent experiments [37,38].Here we combine θ r radial profiles of cases {I}-{IV} with corresponding 2D mode structures to illustrate their underlying connections.Since the poloidal harmonics of e-BAE are weakly coupled, we choose the dominant m = 6 harmonic of δφ and calculate θ r according to δφ m=6 = |δφ m=6 |exp(iθ r ).In figure 8, the blue solid line represents θ r , the red solid and red dashed lines represent the real and imaginary parts of δφ m=6 respectively, and the radial domain of e-BAE eigenfunction with finite amplitude is marked as a gray shaded region.As shown in figures 8 (a) and (c), θ r profiles almost remain constant in the e-BAE region of cases {I} and {III}, which correspond to the relatively symmetric mode structures in figures 7 (a1) and (c1) respectively.In contrast, θ r profiles show large variations in the e-BAE region of cases {II} and {IV} due to the anti-Hermitian EE-KPC term that breaks the poloidal up-down symmetry.However, the poloidal phase shifts at different radial locations are approximately symmetric with respect to the e-BAE peak as shown in figures 8 (b) and (d) (because the EE gradient profile peaks at the mode rational surface in our simulations), which still indicate the radial symmetry of 'boomerang' shape mode structures characterized with two symmetric tails as shown in figures 7 (b1) and (d1).In addition, we show the real frequencies, radial positions and mode widths of cases {I}-{IV} on corresponding n = 3 continuous spectra in figure 9.   To further investigate e-BAE frequency and growth rate dependencies on EE-drive, we vary n h0 and T h0 in amplitudes while remaining profiles unchanged.The e-BAE frequency exhibits weak dependence on n h0 and T h0 , which slowly decreases as n h0 and T h0 increase in figure 10.In contrast, the growth rate is sensitive to both n h0 and T h0 , which linearly increases with n h0 due to the enhanced drive in figure 10 (a), and initially increases and then decreases with T h0 in figure 10 (b) according to the requirement of resonance condition.The intersections of the black dashed lines and growth rate curves in figure 10 represent the marginal stable e-BAEs, of which n h0 and T h0 values give the excitation thresholds that EE-drive just overcomes the continuum damping [39,40] and radiative damping/Landau damping of bulk plasma [41,42,43].From figures 6 and 10, it is seen that both the mode structure and dispersion relation of e-BAE from MAS simulations quantitatively agree with the first-principle GTC results in figures 2 and 3 of Ref. [16], where the deeply-trapped EEs play a dominant role for e-BAE excitation.
Moreover, in order to clarify the competition between EE excitation and bulk plasma damping, we compute the EE drive γ drive in the absence of dissipation by droping Landau damping terms associated to |k || | in Eqs. ( 25), ( 26) and ( 27), and compute the bulk plasma damping rate γ damp by isolating EE drive effect (i.e., only keeping the real part of EE response functions in Eqs. ( 21) and ( 22) while droping the imaginary part), and γ drive and γ damp dependencies on n h0 and T h0 are shown in figure 11.It should be noted that both bulk plasma and EE non-perturbative effects, namely, modifications on mode structure and real frequency, are kept for calculating γ drive and/or γ damp .For example, γ damp shows relatively strong dependency on T h0 in figure 11 (b), which indicates that EE non-perturbative effect can quantitatively affect the bulk plasma damping on e-BAE through altering the mode structure and real frequency, and thus becomes necessary for accurate damping and growth rate calculation.Meanwhile, the γ drive − |γ damp | agrees well with the gross growth rate γ using comprehensive model with both EE drive and bulk plasma damping, which confirms the correctness of our method for obtaining γ drive and γ damp that covers not only the unstable regime of |γ damp | γ drive but also the marginally stable regime of |γ damp | ∼ γ drive , and thus provides useful damping and drive information for linear stability and further nonlinear dynamics analyses [44].

Conclusions
In summary, we have formulated a novel fluid-kinetic hybrid model that couples drift-kinetic EEs to the Landau-fluid bulk plasmas in general geometry, which retains the deeply-trapped EE physics in the simulations of kinetic-MHD processes in a non-perturbative manner.The new model has been implemented and verified in eigenvalue code MAS [24] with practical applications that cover MHD modes, AEs and drift wave instabilities, and the main characteristics are summarized as follows.
1. Drift-kinetic description of EE responses.The EE perturbed distribution is solved from the driftkinetic equation with well-circulating approximation for passing EEs and deeply-trapped approximation for trapped EEs, which keeps the EE kinetic effect of precession drift resonance that is responsible for the excitations of most EE-driven AEs, and the EE fluid effects such as adiabatic and convective responses.Although the well-circulating and deeply-trapped approximations are made in the derivation, it is shown that the EE moments integrated from the perturbed distribution, i.e., perturbed density, parallel current and pressure, can well guarantee the conservation property of EE continuity equation.

2.
Improved deeply-trapped model.The deeply-trapped approximation is applied for deriving the EE-drive terms with dominant precessional drift resonance, which has computational advantage because the heavily numerical integration along the realistic particle orbit is not involved.Specifically, this approximation is made for both the calculation of precession frequency (i.e., ω d ) and bounce-average operations on electromagnetic fields (i.e., δφ, δψ and ω d δψ), where we induce a control parameter λ low to calculate the 2D poloidal profile of deeply-trapped particle fraction f t , namely, only keep the trapped particles that satisfy ω d /ω d0 ∼ 1 and |θ| 1/|nq − m| simultaneously.By improving the accuracy of deeply-trapped approximation with above two constraints, MAS simulations of e-BAE mode structure and dispersion relation show good agreements with gyrokinetic PIC simulations [16].
3. Non-perturbative approach.The EEs are self-consistently incorporated in MAS computations of mode structure, real frequency and growth rate.The non-perturbative effects of EEs on e-BAE mode structures and corresponding radial variations of phase angle profiles are demonstrated.In particular, the EE-KPC term is found to effectively twist the e-BAE global mode structure and break the poloidal up-down symmetry due to the anti-Hermitian contributions to the dielectric tensor, while EE-IC term represents the fluid-like convective response and thus compensate the Hermitian part, which leads to the mode structure to be more symmetric.The EE-KPC term is also responsible for the large radial variations of phase angle, of which poloidal phase shifts at different radial locations are consistent with the 'boomerang' shape mode structure.

4.
Comprehensive damping and driving effects.The original Landau-fluid model for bulk plasmas in MAS [24] has already incorporated the important continuum damping, Landau damping and radiative damping.The newly formulated fluid-kinetic hybrid model combines the EE-drive together with various damping mechanisms, which is more reliable for explaining the marginal stable AEs observed in experiments and is useful for evaluating the EE density and temperature thresholds for AE excitation.
For comparison, the bulk plasma dampings are commonly absent in traditional kinetic-MHD hybrid simulations that might affect the unstable spectra.
With the high efficiency of eigenvalue approach in computation, the upgraded MAS is suitable for fast parameter scans of deeply-trapped EE relevant problems that attract attentions in experiments.Meanwhile, the coupling scheme of EE and bulk plasma can be used for other distributions and/or numerical integration along realistic orbits in phase space, e.g., showing-down distribution and barely-trapped EEs.
ω d ω b for trapped EE particles, where ω b = 2π/τ b and τ b = dl/v || is the bounce period, and ω d = ω d dl/v || /τ b is the precession frequency, and l is the traveled distance of trapped particle along magnetic field line in one bounce motion period.Defining ω b ∂/∂η = v || b 0 • ∇ for trapped particle with η being the bounce angle coordinate and dη = ω b dl/v || , we can transform Eq. (6) into banana orbit center frame using δK t = δK t b exp (iα) with α = − dη (ω d − ω d ) /ω b

)Figure 1 :
Figure 1: The schematic diagram of fluid-kinetic hybrid model.The solid boxes represent the formulation used in computation, while the equations in the dashed boxes are only used for the model derivation.The arrow lines indicate that the vorticity equation is derived by combining the thermal ion, thermal electron and EE continuity equations through the quasi-neutrality condition and parallel Ampere's law.The red circles label the unknown physical variables calculated in corresponding coupling equations.
) for comparison, while δn N A h and δP N A h in Eqs.(21) and (22) contain the response functions R n √ ζ and ζR n √ ζ as shown in figure 2, and we use the limit values of δn N A h and δP N A h at ζ → 0 and ζ → +∞ to compare with thermal electrons, which read

Figure 3 :
Figure 3: The safety factor q profile and magnetic shear s = 1 q

Figure 4 :
Figure 4: Analytic geometry using concentric circular.(a) The ratio between the exact precession frequency ω d (ψ b , θ b ) and the deeply-trapped approximation ω d0 (ψ b ) of all trapped EEs in the poloidal plane.Note that ω d /ω d0 is plotted according to banana tip coordinates (ψ b , θ b ).The trapped particle fraction f t in Eq. (19) with different lower cutoffs of pitch angle (b) λ low = B a /B max and (c) λ low = 1.

Figure 5 :
Figure 5: DIII-D general geometry.The captions of (a)-(c) are the same with figure 4.

Table 1 :
Four simulation cases of n = 3 e-BAE with different EE physics.Note that EE-IC and EE-KPC terms are given in Eq. (24), which represent the fluid and kinetic EE responses respectively.Case EE-IC EE-KPC ω r (V Ap /R 0 ) γ (V Ap /R 0 )

Figure 7 :
Figure 7: The 2D poloidal mode structures of electrostatic potential δφ with different EE terms in Eq. (24).(a) Case (I): drop both EE-IC and EE-KPC terms.(b) Case (II): drop EE-IC term and keep EE-KPC term.(c) Case (III): keep EE-IC term and drop EE-KPC term.(d) Case (IV): keep both EE-IC and EE-KPC terms.

Figure 8 :
Figure 8: (a)-(d) The dominant m = 6 harmonic radial profiles of electrostatic potential δφ and corresponding phase angle θ r for cases (I)-(IV).The radial profile of θ r = arctan(Im(δφ m=6 )/Re(δφ m=6 )) is indicated by the blue solid line, and the radial profiles of Re (δφ m=6 ) and Im (δφ m=6 ) are indicated by the red solid line and red dashed line respectively.The gray shaded region indicates the radial domain of δφ m=6 with high intensity.

Figure 9 :
Figure 9: The e-BAE frequencies of cases (I)-(IV) and continuous spectra.The thick lines represent the Alfvénic branch and the thin lines represent the acoustic branch.

Figure 10 :
Figure 10: The comparisons of e-BAE frequency and growth rate between MAS (blue and red colors) and GTC (black color) simulations.(a) Scan the EE density n h0 with fixed T h0 = 25T e0 .(b) Scan the EE temperature T h0 with fixed n h0 = 0.05n e0 .

Figure 11 :
Figure 11: The kinetic EE drive γ drive (the magenta squares) and bulk plasma damping γ damp (the green squares) calculated by MAS.(a) Scan the EE density n h0 with fixed T h0 = 25T e0 .(b) Scan the EE temperature T h0 with fixed n h0 = 0.05n e0 .The individual calculations of γ drive and γ damp are consistent with the gross growth rate γ (the red squares) obtained from simulations that incorporate both EE drive and bulk plasma damping.