Extended hybrid model

# An extended hybrid magnetohydrodynamics gyrokinetic model for numerical simulation of shear Alfvén waves in burning plasmas

## Abstract

Adopting the theoretical framework for the generalized fishbonelike dispersion relation, an extended hybrid magnetohydrodynamics gyrokinetic simulation model has been derived analytically by taking into account both thermal ion compressibility and diamagnetic effects in addition to energetic particle kinetic behaviors. The extended model has been used for implementing an eXtended version of Hybrid Magnetohydrodynamics Gyrokinetic Code (XHMGC) to study thermal ion kinetic effects on Alfvénic modes driven by energetic particles, such as kinetic beta induced Alfvén eigenmodes in tokamak fusion plasmas.

## 1 Introduction and motivation

Nonlinear numerical simulations of magnetohydrodynamics (MHD) and Alfvén modes driven by energetic particles (EPs) mostly rely on hybrid MHD gyrokinetic codes, such as HMGC [1], M3D [2], and MEGA [3]. In the hybrid MHD gyrokinetic model, the thermal plasma component is described by MHD, while EP dynamics, in the so-called pressure coupling equation [2], is accounted for via the divergence of the EP pressure tensor, which is computed by solving the gyrokinetic equation with particle in cell (PIC) techniques. Kinetic treatments of the thermal plasma component are well known and generally implemented in (linear) spectral codes, such as NOVA-K [5, 6] and MARS-K [7]. More recently, significant developments in gyrokinetic simulation codes, such as GTC [8] and GYRO [9], have also allowed investigating the kinetic effects of thermal plasma and EP dynamics on long wavelength electro-magnetic fluctuations, which were previously investigated only with codes based on the hybrid MHD-gyrokinetic approach [1, 2, 3, 4]. In HMGC [1], the thermal plasma description is originally limited to the reduced MHD model [10]. In the present work, our goal is to extend the hybrid MHD-gyrokinetic model implemented in HMGC [1] to the low-frequency domain of the beta induced Alfvén eigenmode (BAE) - shear Alfvén wave (SAW) continuous spectrum [11], where the mode frequency can be generally comparable with thermal ion diamagnetic and/or transit frequencies, i.e. . In this frequency range, where kinetic thermal ion (KTI) gap generally exists and influences plasma dynamics [11], there is a continuous transition between various MHD and SAW fluctuation branches, as predicted theoretically [12, 13, 14, 15, 16, 17, 18, 19, 20] and confirmed experimentally [21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. Another notable feature of these low frequency fluctuations is that they may be resonantly excited by wave-particle interactions with EPs as well as thermal plasma particles, depending on the perpendicular wavelength [20, 31]. With the extended hybrid MHD gyrokinetic model discussed here, it will be possible to investigate various problems related with resonant excitation of Alfvénic and MHD fluctuations by EPs in the BAE-SAW continuous spectrum, consistent with gyrokinetic codes, e.g., GTC [32], in a common validity domain. Therefore, both the eXtended HMGC (XHMGC) and GTC codes can be verified using different models; yielding more detailed understanding of the underlying physics. In fact, theoretical and numerical work, presented in this article and partly developed within the framework of the SciDAC project on “Gyrokinetic Simulation of Energetic Particle Turbulence and Transport” (GSEP), was the prerequisite for successful verification of XHMGC predictions against analytic theories [33] as well as GTC numerical simulation results [34, 35], reported recently.

In this work, we extend the hybrid MHD-gyrokinetic model, derived originally in [2] for applications to numerical simulations of EP driven Alfvén modes. The main differences with respect to the usual pressure coupling equation [2] are due to renormalization of the inertia term, to properly account for finite thermal ion diamagnetic effects, as well as to the gyrokinetic treatment of the thermal ion pressure tensor, which allows us to properly handle wave-particle resonant interactions in the low frequency regime, where they can be of crucial importance for the analysis of linear and nonlinear behaviors of collisionless burning plasmas. The extended model has been developed assuming ideal Ohm’s law as well as ignoring finite Larmor radius (FLR) effects in order to simplify the technical complications while still maintaining all essential physics ingredients [36]. In practice, maintaining the ideal MHD Ohm’s law as limiting case implies assuming and neglecting ion FLR effects, although finite magnetic drift orbit widths (FOW) are fully retained [37]. A more general approach without these simplifying assumptions will be developed in a separate work. For demonstrating the validity of the modified equations, we show that they are equivalent to the quasi-neutrality and vorticity equations derived in [36] for the frequency range from the kinetic ballooning mode (KBM) and BAE to the toroidal Alfvén eigenmode (TAE). The XHMGC model equations in the linear limit are equivalent to the extended kinetic MHD used in spectral codes, such as NOVA-K [5, 6] and MARS-K [7], but with EP dynamics treated non-perturbatively and on the same footing as the thermal plasma response (see Section 2 for more details). The possibility of investigating nonlinear dynamics, however, makes XHMGC more suitable to direct comparisons with M3D [2] or gyrokinetic codes [32, 38] in a common validity domain.

The paper is organized as follows. In Section 2, the extended hybrid model equations are presented and discussed within the theoretical framework of [36]. In Section 3, we describe the numerical implementation of the extended model into HMGC, by adding both thermal ion compressibility and diamagnetic effects (of thermal ions as well as EPs) into MHD equations and a thermal ion population in the PIC module. In Section 4, possible applications and validity limits of XHMGC are discussed. A synthetic summary of current BAE numerical simulation results [33] are also provided. Finally, conclusions and discussions are given in Section 5.

## 2 Derivation of the extended hybrid model

Reference [36] presents a general theoretical framework for stability analyses of various modes and the respective governing equations. It shows that all modes of the shear Alfvén branch having frequencies in the range between the thermal ion transit and Alfvén frequency can be consistently described by one single general fishbone-like dispersion relation (GFLDR) [12, 16, 36, 15, 14]. Reference [36] discusses various reduced equations governing the evolution of SAW fluctuations in burning plasmas, using the general approach of reference [39]. In this sense, reference [36] could seem not the optimal eliminate reference framework for further generalizing the HMGC hybrid model equations  [1, 40], which are to be used for nonlinear studies as well. However, the detailed analyses of reduced model equations, reported in [36], on the basis of specializations of ordering of dimensionless parameters in the case of burning plasmas of fusion interest, starting from the somewhat different orderings of interest to space plasmas given in [39], allow to fully grasp the physics implications of the underlying approximations. Moreover, on the basis of our discussions, it is straightforward to motivate the extension of the derived model equations to the nonlinear case, as shown at the end of this section.

Considering that the characteristic frequency, , is much lower than the ion cyclotron frequency, , we may adopt the gyrokinetic theoretical approach and closely follow reference [39]. The low-frequency plasma oscillations can, thus, be described in terms of three fluctuating scalar fields: the scalar potential perturbation , the parallel (to , with the equilibrium magnetic field) magnetic field perturbation and the perturbed field , which is related to the parallel vector potential fluctuation by

 δA∥≡−i(cω)b⋅∇δψ. (1)

The governing equations for describing the excitation of the shear Alfvén frequency spectrum by energetic ions precession, precession-bounce and transit resonances in the range , covering the entire frequency range from KBM/BAE [14, 13, 21, 41] to TAE [42, 43, 44], are generalized kinetic vorticity equation and quasi-neutrality condition, which can be written as followings, in the limit of vanishing FLR (see equation (16) and equation (17) in reference [36]):

 B0⋅∇(k2⊥k2θB20B0⋅∇δψ)+ω(ω−ω∗pi−nEmEnimiω∗pE)v2Ak2⊥k2θδϕ −⟨∑s≠e4πesk2θc2ω^ωdsδKs⟩+∑s4πk2θB20k×b⋅∇(Ps⊥+Ps∥)Ωκδψ=0, (2)
 ⟨∑s≠Ee2sms∂F0s∂ε⟩(δϕ−δψ)+∑s=ies⟨δKs⟩=0, (3)

where the non-adiabatic particle response, , is obtained via the drift-kinetic equation

 [ωtr∂θ−i(ω−ωd)]sδKs=i(em)sQF0s[(δϕ−δψ)+(^ωdω)sδψ]. (4)

Here, angular brackets stand for velocity space integration, denotes all particle species ( bulk electrons, bulk ions, energetic particles), and are the species electric charge and mass, is the equilibrium distribution function (generally anisotropic), the energy per unit mass, , , is the wave vector, is the cyclotron frequency, is the perpendicular wave vector, is the diamagnetic frequency, and are, respectively, the total perpendicular and parallel plasma pressures, is the transit frequency and , with and . Note that the difference between and , with , has been discussed in [36, 39] and, generally, must be handled properly; although, for many applications in low pressure () plasmas, one can consider after solving for from perpendicular pressure balance [39, 36], as implicitly assumed in equations 23 and 4. Note, also, that we have maintained the EP contribution to the divergence of the polarization current, which is represented by its leading term in equation 2. This term is readily derived from the last term on the left hand side (LHS) of equation (13) in reference [36] (see also A for further details) and was neglected in there due to the ordering , valid in a burning plasma dominated by fusion alpha particle self-heating. Here, and denote the beta values of EP and bulk plasma components (electrons and thermal ions), respectively, while and are the energy confinement time and EP slowing down time. More generally [12, 13, 14, 15], the ordering better represents nowadays magnetized plasmas of fusion interest and, thus, , as assumed in equation 2.

Equations 2 and 3, together with the drift-kinetic equation, equation 4, are the simplest yet relevant equations for analyzing the resonant excitations of SAW by EPs. Equation 2 demonstrates that both resonant as well as non-resonant responses due to the term enter via the magnetic curvature drift coupling. In the high frequency case, , the thermal ion kinetic compression response can be neglected. Thus, the quasi-neutrality condition, equation 3, reduces to the ideal MHD approximation, ; i.e.  [39]. Meanwhile, neglecting the terms, equation 2 becomes equivalent to equation (3) in [2], i.e. the following pressure coupling equation in the hybrid MHD-gyrokinetic approach

 ρbdvbdt=−∇Pb−(∇⋅PE)⊥+J×Bc; (5)

where the subscript denotes the bulk plasma (electrons and thermal ions), while and are, respectively, bulk plasma mass density and fluid velocity. Here, the EP contribution to the perpendicular momentum change of the plasma has been neglected, due to  [36, 2], and thermal ion diamagnetic effects are consistently dropped since .

In order to extend the hybrid model to the low-frequency regime where , we need to include the effects of the thermal ion compressibility within the hybrid simulation scheme. That is, we need to include effects associated with the terms in equations 2 and 3. First, in order to simplify the discussions, we formally assume in the present work; the general case with finite will be considered elsewhere. Thus, according to equation 3, we have and the ideal MHD condition remains valid. Next, we proceed to establish correspondences between the pressure coupling equation, equation 5, and the generalized kinetic vorticity equation, equation 2.

Applying the operator to the linearized equation 5 and noting the quasi-neutrality condition , we readily derive

 1c∂∂tB0⋅∇δJ∥B0i+1c∂∂tδB⋅∇(J∥0B0)ii+∂∂t∇⋅⎛⎜⎝B0×ρb0dδvbdtB20⎞⎟⎠iii +∂∂t∇δPb⋅(∇×bB0)iv+∂∂t∇⋅(b×(∇⋅δPE)⊥B0)v=0. (6)

Noting also the parallel Ampère’s law along with ,

 4πδJ∥=−c∇2δA∥, (7)

and equation 1, term (i) can be seen to correspond to the field line bending term; i.e. the first term in equation 2. Term (ii), on the contrary, does not have any direct correspondence in equation 2. This term is the usual kink drive and it was dropped in the analysis of [36], focusing on drift Alfvén fluctuations with high mode numbers, for it is formally of , with the toroidal mode number. However, as noted in equation (A1) of [36], term (ii) is readily recovered in a form that can be straightforwardly reduced to that reported here. Meanwhile, from the linearized Ohm’s law

 δE⊥+1cδvb×B0=0, (8)

and , term (iii) corresponds to the second term in equation 2 with the terms neglected. To establish correspondences between the pressure responses in equations 2 and 2, we first denote . It can then be shown (see appendix B) that term (iv) corresponds to the thermal ion and electron contributions to the last term on the LHS of equation 2, when kinetic compression effects of the background thermal plasma are neglected.

Finally, let us discuss term (v), due to EP pressure perturbation, which can be expressed as (see C).

 ∂∂t∇⋅(b×(∇⋅δPE)⊥B0)=ωB0Ωκ(δPE∥+δPE⊥). (9)

Meanwhile, noting the definition of  [39, 36], the term in equation 2 can be shown to be related with the pressure perturbations as (see D)

 ⟨4πesk2θc2ω^ωdsδKs⟩ = 4πωk2θcB0Ωκ(δPs⊥+δPs∥) (10) +4πk2θB20(k×b)⋅(∇P0s⊥+∇P0s∥)Ωκδψ.

Note that the 2nd term in the right hand side (RHS) disappears in the ideal MHD limit. Equation 10, thus, clearly demonstrate that the contribution in equation 2, combined with the 3rd term on the RHS of equation 10 (or the last term on the LHS in equation 2), has the same form of equation 9 and recovers the total pressure response of term (iv) in equation 2 for . In other words, the term in equation 2 corresponds to the kinetic compressibility component of the pressure perturbations.

Summarizing the above discussions, it is clear that, in order to include effects due to finite thermal ion compressibility and diamagnetic drift as well as the finite EP contribution to the divergence of the polarization current, the pressure coupling equation in the MHD-gyrokinetic approach, equation 5, has to be modified such that its perpendicular components are given by equation 26 of A, which we rewrite here for the reader’s convenience

 [ρb(∂∂t+vb⋅∇)+b×∇P0E⊥ωcE⋅∇]δvb= −∇⊥Pe−(∇⋅Pi)⊥−(∇⋅PE)⊥+(J×Bc)⊥. (11)

Here, , and the “unshifted” pressure tensors and need to be calculated from solutions of the gyrokinetic equations as specified in A, while is consistently neglected in the present approach, assuming . Reminding the concluding remark of A, this equation readily reduces to the well-known pressure coupling equation 2, in the limit where thermal ion diamagnetic effects and EP contribution to the divergence of the polarization current are neglected.

As anticipated above, in the present work, we followed reference [36], since that has a detailed discussion of validity limits of different reduced models of the whole vorticity and quasi-neutrality equations, derived for fusion applications and following the trace of reference [39]. Equation 2 includes equilibrium parallel current effects, as discussed earlier in this section and in [36] (Appendix). This simple remark readily follows from the discussion presented in [34] as well as the modified momentum balance equation implemented in XHMGC, i.e. equation 2 itself. The present model is valid in the nonlinear case too, as shown by the simple derivation provided in A and by the following discussion. This is deduced easily from direct inspection of equation (5) in reference [45]. That equation clearly shows that, for the small FLR limit considered in HMGC [1, 40], the nonlinear terms, treated explicitly, are those that are coming from convective nonlinearity and from the Maxwell stress nonlinearity, when the thermal ion response is taken in the fluid limit, both of which are readily obtained from equation 2 upon application of the operator , as it was done for equation 5 earlier in the section. Other nonlinear dynamics, which are implicitly included in and terms, are fully retained via equation 2. Thus, the back reaction of zonal structures onto SAW fluctuations is fully accounted for, i.e. that of zonal flows (ZFs) and fields as well as radial modulation of equilibrium profiles [9, 11, 46] which also enter via the diamagnetic terms in equation 2, computed on the whole (slowly evolving) thermal ion and EP pressure profile, obtained from the respective toroidally and poloidally averaged distribution functions. This choice is consistent with known approaches to nonlinear MHD equations, accounting for finite diamagnetic drift corrections [47, 48, 49, 50, 51].

Thus, the approximations involved with the extended implementation within XHMGC on the basis of equation 2 consist of neglecting FLR, assuming electron as a massless fluid, considering (such that parallel Ohm’s law is recovered in the ideal MHD limit) and accounting for Reynolds stress in the thermal ion fluid limit. The possible further extension of the present model to include finite and generalizing the parallel Ohm’s law, while maintaining other simplifying assumptions, is straightforward on the basis of the present discussion and will be reported in a separate work [52]. Here we note that the present extended hybrid model, based on equation 2, with clearly formulated assumptions that limit its applicability, includes very rich physics; e.g. it is capable to correctly evaluate the renormalized inertia for ZFs, for which the trapped thermal ion dynamics is of crucial importance, and to account for geodesic acoustic mode (GAM) kinetic response, including Landau damping.

So far, XHMGC has been used for moderate EP drive [33, 34, 35, 37], where the EP diamagnetic correction to the divergence of the polarization current can be neglected, as argued in [36]. Actually, in the studies reported in [33], thermal ion diamagnetic contribution to the polarization current is also neglected, since the case of uniform thermal ion pressure profiles is investigated in there for facilitating comparisons of numerical simulation results with analytic theory predictions (see also section 4).

## 3 Numerical implementation

HMGC [1] is used for investigating linear and nonlinear properties of moderate toroidal number (n) shear Alfvén modes in tokamaks. It solves the coupled set of O() reduced-MHD equations [10] for the electromagnetic fields and the gyro-center Vlasov equation for a population of energetic ions, where large aspect ratio is assumed, i.e. , with and the tokamak minor and major radius, respectively . Energetic particles contribute to the dynamic evolution of the wave fields via the pressure tensor term in the MHD equations, as described by the pressure coupling equation [2]. This code allows us to describe both self-consistent mode structures in toroidal equilibria and EP dynamics, as well as to get a deeper insight into how the Alfvénic modes affect the confinement of such particles.

The extended model, described in section 2, has been implemented into the eXtended version of HMGC (XHMGC). Following the general procedure, described in references [1, 40], for the formal manipulation of equation 2, the relevant equations for the MHD solver are in terms of the poloidal magnetic field stream function and , which is proportional to the scalar potential and defined as , can be written in the following form in the cylindrical coordinate system :

 ∂Ψ∂t=R2R0∇Ψ×∇φ⋅∇U+B0R0∂U∂φ+ηc24π△∗Ψ+O(ϵ4vABφ), (12)
 ^ρ(DDt+2R0∂U∂Z)∇2⊥U+∇^ρ⋅(DDt+1R0∂U∂Z)∇⊥U (13) −(R2R30ωci0∂P0i⊥∂Z+R2R30ωcE0∂P0E⊥∂Z)∇2⊥U −∇(R2R30ωci0∂P0i⊥∂Z+R2R30ωcE0∂P0E⊥∂Z)⋅∇⊥U +∇⋅[R4R30(∇φ×∇P0E⊥ωcE0+∇φ×∇P0i⊥ωci0)⋅∇∇⊥U] = 14πB⋅∇△∗Ψ+1R0∇⋅[R2(∇Pe+∇⋅Πi+∇⋅ΠE)×∇φ] +O(ϵ4ρv4Aa2),

where we have maintained the same notation of reference [1] and explicitly show the additional terms that have been added to implement the extended XHMGC model. Thus,

 ^ρ=R2R20ρ,  DDt=∂∂t+R2R0∇U×∇φ⋅∇, ∇2⊥≡1R∂∂RR∂∂R+∂2∂Z2,

the Grad-Shafranov operator is defined by

 △∗≡R∂∂R1R∂∂R+∂2∂Z2, (14)

is the vacuum magnetic field on the magnetic axis at  1, , and the subscript denotes components perpendicular to . In the above equations, is the fluid velocity, is the bulk plasma mass density, is the scalar pressure of bulk electrons, and are, respectively, the pressure-tensor of the EP and thermal ions, computed with the definitions of equations 20 and 24, given in A, and, is the resistivity and is the speed of light. These equations have been first derived in reference [10], limited to the MHD description of the thermal plasma, while the inclusion of energetic particle dynamics has been discussed in references. [1, 40]. Here, in the proposed further extension of the numerical model, thermal ion dynamics as well as diamagnetic effects are also taken into account, according to equation 2, derived in A and section 2. At the leading order in , , the reduced-MHD equations describe the thermal plasma in the cylindrical approximation. Toroidal geometry enters the equations as corrections at the next order in the inverse aspect ratio.

In order to close equations 12 and 13, the EP and thermal ion pressure tensor components can be obtained by directly calculating the appropriate velocity moments of the distribution function for the particle population interacting with the perturbed electromagnetic field. As discussed in section 2, we initially assume the limit for the sake of simplicity, i.e. . Meanwhile, with cold electron assumption and ignoring thermal ion finite Larmor radius (FLR), ideal MHD parallel Ohm’s law can be readily recovered.

As to numerical formulation, the equations of motion in gyro-center coordinates for thermal ions are in the same form, mutatis mutandis, as those reported in [1] for EPs. In the gyrocenter-coordinate system , where is the gyrocenter position, is the conserved magnetic moment, is the parallel speed and is the gyrophase, the equations of motion take the form

 d¯Rdt = ¯Vb+esmsΩsb×∇ϕ−¯VmsΩsb×∇a∥ +[¯Mms+¯VΩs(¯V+a∥ms)]b×∇lnB, d¯Mdt = 0, d¯Vdt = 1msb⋅{[esΩs(¯V+a∥ms)∇ϕ+¯Mms∇a∥]×∇lnB (15) +esmsΩs∇a∥×∇ϕ}−ΩE¯Mmsb⋅∇lnB.

Here, the subscript denotes either EP or thermal ion species and, using the same notations as in [1], is the corresponding cyclotron frequency. The fluctuating potential is related to the poloidal magnetic field stream function through the relationship . The parallel electric field term in the equation for has been suppressed, neglecting, thus, small resistive corrections to the ideal-MHD parallel Ohm’s law. Meanwhile, the pressure tensor can be written, in terms of the gyrocenter coordinates, as

 Πs(t,x)=1m2s∫d¯ZDZc→¯Z¯Fs(t,¯R,¯M,¯V) [Ωs¯MmsI+bb(¯V2−Ωs¯Mms)]δ(x−¯R), (16)

where is the unit tensor, , is the gyrocenter distribution function and is the Jacobian of the transformation from canonical to gyrocenter coordinates. The distribution function satisfies the Vlasov equation

 (∂∂t+d¯Rdt⋅∇+d¯Vdt∂∂¯V)¯Fs=0, (17)

where and are given by equation 3. In the numerical implimentation of XHMGC, equations 3 and 17 can be readily solved as a full-F simulation. On the other hand, a algorithm [40, 53, 54, 55] is also implemented in order to minimize the discrete particle noise. The latter is recommended as far as , the former when .

## 4 Applications

In general, the extended version of HMGC can have two species of kinetic particles. On one hand, one can use XHMGC for investigating thermal ion kinetic effects on Alfvénic modes driven by EP. On the other hand, it may be interesting to use XHMGC as a tool to simulate two coexisting EP species, generated e.g. by both ion cyclotron resonance heating (ICRH) and neutral beam injection (NBI) heating, in order to study linear excitation of Alfvénic fluctuations and Energetic Particle Modes (EPM) [15], as well as the interplay between the respective nonlinear physics controlled by the different heating sources [56].

HMGC has been extensively used in [1] and [57] to investigate the linear physics (damping and EP drive mechnisms), and in [40] and [58] to analyze the nonlinear dynamics of EPM. XHMGC has been verified against those previous findings and can recover numerical simulation results in the above studies. Furthermore, by accounting for the kinetic thermal ion effects, XHMGC shows the existence of Kinetic BAE (KBAE) which can be seen as radially trapped eigenstates due to discretization of BAE-SAW continuum by FLR/FOW effects, as well as KBAE resonantly excited by wave-particle interactions with EPs [33, 37].

As an example to demonstrate the capability of XHMGC, we briefly report simulation results of KBAE, which is discussed in detail in [33]. The results show that a fully kinetic treatment of thermal ions is necessary for a proper description of the low frequency Alfvénic fluctuation spectrum. By including thermal ion compressibility, our numerical simulations do show the existence of a finite-frequency BAE accumulation point in the SAW continuum, which was demonstrated analytically and numerically using MHD codes [59]. Meanwhile, when effects due to finite ion drift orbit width (FOW) are included, our simulations clearly demonstrate that the BAE-SAW continuum becomes discretized; yielding a series of discrete kinetic eigenmodes with small frequency separation [33, 60]. In figure 1, we have plotted the BAE accumulation frequencies in the fluid limit which are defined as with in the current case, as well as eigenmode frequencies obtained from simulation results. The analytically predicted KBAE frequencies [33] are in good agreement with observations from numerical simulations. The results also indicate that FOW kinetic effects increase with the toroidal mode number, as expected [33, 60].

On the other hand, our simulations also show that KBAE can be driven by EPs. In figure 2, we can see that the frequencies scale properly with the KBAE frequencies; and the growth rates decrease with the thermal ion temperature due to the stronger ion Landau damping and/or the weaker EP drive due to the increased frequency mismatch between mode and characteristic EP frequencies. In the absence of thermal ion kinetic effects, the excited modes may be identified as energetic particle mode (EPM); which requires sufficiently strong drive to overcome the SAW continuum damping. Including the thermal ion kinetic effects not only introduce a finite kinetic thermal ion frequency gap at the BAE accumulation frequency but also discretize the BAE-SAW continuum. In that case, the continuum damping is greatly reduced or nullified, and the discrete KBAE’s are more readily excited by the EP drive.

## 5 Conclusions and discussions

In the present work, we have employed the theoretical framework (generalized kinetic vorticity and quasi-neutrality equations) of the generalized linear fishbone dispersion relation and derived an extended hybrid MHD-gyrokinetic simulation model applicable to the low-frequency regime, where effects of thermal ion compressibility and diamagnetic drifts play significant roles in the dynamics of Alfvén waves and energetic particles in tokamak plasmas. The extended simulation model has been implemented into an eXtended version of HMGC (XHMGC). Initial simulations of XHMGC have discovered the existence of KBAE discretized by the thermal ion FOW effects, which are absent in conventional MHD codes. Simulations also demonstrate that KBAE can be readily excited by EPs. In the current model, we have taken and neglected finite Larmor radius effects in order to simplify the presentation and focus on the most important qualitative new physics connected with implementation of the thermal ion compressibility. In addition, XHMGC is limited to circular shifted magnetic surfaces equilibria, with relatively large aspect ratio; XHMGC includes kinetic effects related to both bulk and fast ions; however, it is typically used for retaining only the perturbed pressure for two EP species; XHMGC doesn’t include rotation (see A), while it retains the perturbed electrostatic potential. These additional effects will be considered in future works.

More recently, the electromagnetic formulation [8] of global gyrokinetic particle simulation in toroidal geometry has been implemented in GTC [32]. In such a code, ions are treated by the gyrokinetic equation, while electrons are simulated using an improved fluid-kinetic electron model [8]. In [34], the connection between the extended hybrid MHD-gyrokinetic model and gyrokinetic simulation model has been verified in the drift kinetic limit as well as ignoring the terms on the order of . Instead of directly calculating the pressure tensor, lower moments of the kinetic equation have been calculated, i.e. the perturbed density and parallel current. Using charge neutrality condition, it can be demonstrated that the combination of the perturbed density and parallel current contribution is totally equivalent to the pressure tensor in equation 2. Therefore, both GTC and XHMGC can be verified using different models in a common validity regime; yielding more detailed understanding of the underlying physics.

On the other hand, kinetic compressibility is also included in (linear) spectral codes, such as NOVA-K [5, 6] and MARS-K [7]. While in NOVA-K only perturbative treatment of the kinetic effects is considered and continuum damping is not taken into account, as in MARS-K, since both are eigenvalue codes, the spectral approach allows the study of the linear stability of all eigenmodes in a general equilibria; meanwhile, kinetic effects are generally related to bulk plasmas only, although the inclusion of fast ions is quite straightforward. As to other hybrid MHD gyrokinetic codes, M3D [2] is based on the pressure coupling equation; MEGA [3, 4] uses a hybrid model for MHD and energetic particles, where the effect of the energetic ions on the MHD fluid is taken into account in the MHD momentum equation through the energetic ion current. The diamagnetic drift effect is evaluated in the MHD equations by adding the diamagnetic advection term to the equation of motion [47, 48, 49, 50, 61]. At present, XHMGC can handle two species kinetic particles self-consistently, but is limited to circular shifted magnetic flux surfaces equilibria with vanishing bulk plasma equilibrium pressure. Meanwhile, a new version of the code with general equilibria is being developed, with the capability of solving fully compressible gyrokinetic particle response.

## Acknowledgments

This work is supported by the ITER-CN under Grant No.2009GB105005, the NSF of China under Grant No. 11075140, Euratom Communities under the contract of Association between EURATOM/ENEA, USDOE GRANTS, SciDAC, and GSEP.

## Appendix A Simple derivation of model equations

Adopting a multi-fluid moment description of plasma dynamics, the force balance equation can be written as

 ρb(∂∂t+vb⋅∇)vb+ρE(∂∂t+vE⋅∇)vE =−∇Pe−∇⋅Pi−∇⋅PE+J×Bc. (18)

Here, and are bulk plasma and EP mass densities, , and from equation (8), having omitted terms that are or higher with respect to the RHS. Furthermore, thermal ion and EP pressure tensors on the RHS have to be interpreted as usual, i.e. with the conventional fluid velocity shift in the definition

 Psij=ms∫dv(vi−usi)(vj−usj)fs, (19)

with the particle distribution function and . When the pressure tensor is computed form the particle distribution function within the gyrokinetic description, some subtleties are connected with the ordering in the plane orthogonal to , with the Larmor radius of the -species, its thermal speed and the characteristic equilibrium radial scale-length. Thus, in the drift-kinetic limit used in this work, , with the unit diagonal tensor and

 Ps⊥ = ms∫dvv2⊥2fs, (20) ^Ps∥ = ms∫dv(v∥−us∥0)2fs. (21)

Note the difference between , used here, and , used in equation 19, being the (slowly evolving) equilibrium particle distribution function. In equation 18, we assumed that only EPs can carry significant parallel fluid velocity.

The perpendicular component of equation 18 can be further simplified, by noting that , with and and the tokamak minor and major radii, and that we are using the optimal ordering . This allows us to rewrite

 [ρE(∂∂t+vE⋅∇)vE]⊥=b×∇P0E⊥ωcE⋅∇vE⊥+ρEu2E∥κ, (22)

where we have dropped the terms, for they are and, similarly, the term, since it is – or, equivalently, – for ; at shorter wavelength or higher frequency, this term would be negligible anyway with respect to the thermal ion inertia response, as negligible would be diamagnetic responses of both EPs and thermal ions. So, equation 22 well describes the physics we want to incorporate in the present analysis. Recalling the definition of , we also have

 −(∇⋅Ps)⊥=−∇⊥Ps⊥−κ(^Ps∥−Ps⊥) (23)

Thus, the second term on the RHS of equation 22 can be combined with the term on the RHS of equation 23, computed for EPs, and actually be reabsorbed into that (up to the relevant order), provided that the pressure tensor is reinterpreted as , with the “unshifted” expression

 Ps∥=ms∫dvv2∥fs. (24)

replacing the usual definition given in equation 21. With this convention on the pressure tensor, the perpendicular components of equation 18 can be rewritten as

 ρb(∂∂t+vb⋅∇)vb+b×∇P0E⊥ωcE⋅∇vE⊥= −∇⊥Pe−(∇⋅Pi)⊥−(∇⋅PE)⊥+(J×Bc)⊥. (25)

Actually, equation 25 can be reduced further when residual terms that are or higher with respect to the RHS are omitted, as noted below equation 18. In fact, one readily obtains

 [ρb(∂∂t+vb⋅∇)+b×∇P0E⊥ωcE⋅∇]δvb= −∇⊥Pe−(∇⋅Pi)⊥−(∇⋅PE)⊥+(J×Bc)⊥. (26)

This equation readily reduces to the well-known pressure coupling equation [2], in the limit where thermal ion diamagnetic effects and EP contribution to the divergence of the polarization current are neglected.

## Appendix B Study of term (iv) in equation 2

In the low- approximation (),

 ∇×(bB0)≅2b×κB0. (27)

Meanwhile, in the incompressible limit,

 ∂∂tδPb+δvb⋅∇P0b=0, (28)

where

 δvb⊥=cB0×∇⊥δϕB20. (29)

Then, with equations 2728 and 29

 ∂t(∇×bB0)⋅∇δPb = 2b×κB0⋅∇∂δPb∂t (30) = c2b×κB0⋅∇(−B0×∇⊥δϕB20⋅∇P0b) = c2b×κB0⋅∇(B0×∇P0bB20⋅∇⊥δϕ) = −2cB20k×b⋅∇P0bΩκδϕ,

where .

## Appendix C Study of term (v) in equation 2

Assuming

 δPE=bbδPE∥+(I−bb)δPE⊥, (31)

we can show

 ∇⋅δPE = (δPE∥−δPE⊥)(b∇⋅b+κ) (32) + b∇∥(δPE∥−δPE⊥)+∇δPE⊥,

where . Thus,

 b×∇⋅PE=b×∇δPE⊥+(δPE∥−δPE⊥)b×κ. (33)

Now

 ∇⋅(bB0×∇δPE⊥)≅2b×κB0⋅∇δPE⊥, (34)

and

 ∇⋅[(δPE∥−δPE⊥)bB0×κ] ≅ b×κB0⋅∇(δPE∥−δPE⊥). (35)

In equation 35, we have used the large aspect ratio assumption, , consistent with the reduced MHD description used in HMGC [1, 40]. Combining equations 33 to 35, we obtain

 ∂∂t∇⋅(b×(∇⋅δPE)⊥B0) = b×κB0⋅∇∂∂t(δPE∥+δPE⊥) (36) = ωB0Ωκ(δPE∥+δPE⊥)

## Appendix D Study of the ∝δKs term in equation 2

Here, we assume the definition of  [36]

 eiLksδKs=δfs−(em)s[∂F0∂εδϕ−QF0ωeiLkJ0(k⊥ρs)δψ]s, (37)

where is the fluctuating particle distribution function, , on the RHS, we have dropped all terms , for they generate contributions of higher order in what follows [36]. Thus, in our treatment, is generally anisotropic, although terms do not appear explicitly. One then finds

 ⟨4πesk2θc2J0(k⊥ρs)ω^ωdsδKs⟩ = ⟨4πesk2θc2ω^ωdsδfs⟩I−⟨4πe2sk2θmsc2ω^ωds∂F0s∂ε⟩δϕII (38) +⟨4πe2sk2θmsc2ω^ωdsQF0sωJ20(k⊥ρs)⟩δψIII,

with denoting velocity integration and is the zero order Bessel function.

For term (I), we obtain

 (I) = 4πωk2θcΩκ⟨ms(μ+v2∥/B)δfs⟩=4πωk2θcBΩκ(δPs⊥+δPs∥), (39)

where

 δPs⊥=⟨m2v2⊥δfs⟩ (40)

and

 δPs∥=⟨mv2∥δfs⟩ (41)

are, respectively, perturbed perpendicular and parallel pressure.

Meanwhile, for , term (III) can be written as

 (III) = ⟨4πesk2θc(μ+v2∥B)Ωκ(ω∂εF0s+mscesB(k×b)⋅