Theory and modeling of electron fishbones

Internal kink instabilities exhibiting fishbone like behavior have been observed in a variety of experiments where a high energy electron population, generated by strong auxiliary heating and/or current drive systems, was present. After briefly reviewing the experimental evidences of energetic electrons driven fishbones, and the main results of linear and nonlinear theory of electron fishbones, the results of global, self-consistent, nonlinear hybrid MHD-Gyrokinetic simulations will be presented. To this purpose, the extended/hybrid MHD-Gyrokinetic code XHMGC will be used. Linear dynamics analysis will enlighten the effect of considering kinetic thermal ion compressibility and diamagnetic response, and kinetic thermal electrons compressibility, in addition to the energetic electron contribution. Nonlinear saturation and energetic electron transport will also be addressed, making extensive use of Hamiltonian mapping techniques, discussing both centrally peaked and off-axis peaked energetic electron profiles. It will be shown that centrally peaked energetic electron profiles are characterized by resonant excitation and nonlinear response of deeply trapped energetic electrons. On the other side, off-axis peaked energetic electron profiles are characterized by resonant excitation and nonlinear response of barely circulating energetic electrons which experience toroidal precession reversal of their motion.


Introduction
The mutual interaction of particle populations, characterized by very disparate kinetic energies, is of great interest for research on thermonuclear plasmas of fusion relevance, and, in particular, for the so-called 'ignited' plasmas, in which the 3.52 MeV α-particles, released in deuterium-tritium (D-T) reactions, have to thermalize by Coulomb collisions with the bulk thermal D-T plasma in order to self sustain its temperature. The interplay of fusion α-particles and magnetohydrodynamics (MHDs), Alfvénic-like modes has been recognized, since long time, as a crucial issue for the success of next generation,'ignited' devices as, e.g., ITER [1]. Indeed, the potential enhancement of the radial transport of energetic particles (EPs) toward the edge of the plasma device while preventing them to fully thermalize could, in turn, degrade the fusion performance on one side, and damage the plasma facing components on the other. Similar phenomenology could also take place because of energetic particles accelerated by auxiliary heating systems, as, e.g., neutral beam (NB) injection and a variety of radio frequency heating and current drive systems, and, indeed, has been observed in a large selection of present days auxiliary heated toroidal plasma devices (see, e.g., [2,3]).
One of the 'case studies' of EP driven MHD-like modes is the 'fishbone' mode, originally observed in the Poloidal Divertor eXperiment (PDX) [4] device, owing its name to the characteristic fishbone-like shape of the perturbed magnetic field signal evolution. The fishbone is an internal kink-like instability driven, in PDX, by energetic ions due to NB injection, which results in anomalous losses of energetic ions themselves. The theoretical interpretation in terms of resonant wave-particle interaction at the EP toroidal precession frequency was first proposed in [5] briefly after the experimental observations. More recently, internal kink instabilities excited by supra-thermal electrons and exhibiting fishbone like behavior have been observed in the Doublet III-D (DIII-D) [6] device, where the high-energy electron population was generated by electron cyclotron (EC) current drive. Later on, other devices observed fishbone oscillations with electron heating only, i.e., electron cyclotron resonant heating (ECRH) and/or lower hybrid heating (LHH) and lower hybrid (LH) current drive (see next section 2 for a brief review of experimental evidence of electron fishbones, e-fishbones). Electron fishbones have linear dispersion relation and excitation mechanisms that are similar to those of energetic ion driven fishbones; moreover, fluctuation induced transport of magnetically trapped resonant particles, due to precession resonance, is expected to depend on energy and not mass of the energetic particles involved, because of the bounce averaged dynamic response [7]. Electron fishbones are characterized by a very small ratio between the resonant particle orbit width and the characteristic fishbone length scale ( x d r , the rigid radial kink-type displacement). This is also expected to be the case of ion fishbones in burning plasmas of fusion interest due to the large plasma current in these devices, while this condition is not realized for the energetic ions in present-day experiments. These analogies between e-fishbones in present-day devices and fishbones in burning plasmas provide a practical motivation for investigating these processes, in addition to the general interest of studying e-fishbones 'per se'.
In this paper, after briefly reviewing the experimental evidences of energetic electrons driven fishbones as observed in present toroidal devices in section 2, and the main results of linear and nonlinear theory of e-fishbones in section 3, the results of global, self-consistent, nonlinear hybrid MHD-Gyrokinetic simulations will be presented in section 4 [8]. In particular, the extended/hybrid MHD-Girokinetic code XHMGC, described in [9], will be used. The effects of considering kinetic thermal ion compressibility and diamagnetic response (in order to allow for an entirely novel treatment of enhanced inertia response [7,10,11] and ion Landau damping [12]), and kinetic thermal electrons compressibility, in addition to the energetic electron contribution, will be enlightened in linear dynamics analysis. Nonlinear saturation and energetic electron transport will also be addressed, making extensive use of Hamiltonian mapping techniques [13,14]. In order to illustrate different nonlinear dynamics, both centrally peaked and off-axis peaked energetic electron profiles will be discussed. In particular, centrally peaked energetic electron profiles are characterized by resonant excitation and nonlinear response of deeply trapped energetic electrons. Meanwhile, barely circulating energetic electrons are identified as responsible of driving the e-fishbone mode and causing its nonlinear evolution for off-axis peaked energetic electron profiles. Final considerations will be given in section 5.

Experimental evidences of energetic electron driven fishbones
Fishbone oscillation driven by energetic ions have been observed for the first time in PDX [4] discharges heated by perpendicular NB injection, where a large n=1 MHD mode was observed causing losses of energetic ions. Deeply trapped ions, in presence of a beam deposition profile peaked near the magnetic axis, were recognized to drive the mode [4,5] because of resonant wave-particle interaction at the EP toroidal precession frequency w d . Fishbone oscillations driven by suprathermal ion population have been observed, since then, on many tokamak devices [2,3,15,16]. Observations indicate that the mode propagates poloidally in the ion diamagnetic drift direction, and toroidally parallel to the EP precession velocity, thus having w w w  r es dh and * * w w w w > ¯0 , consistent with theoretical predictions for unstable modes [5] (see also section 3). Here, w res is the resonance frequency, the overbarx on the quantity 'x' indicates its bounce average, the diamagnetic frequency * * w = =´  · · , with * v s the diamagnetic velocity, v s d the magnetic drift velocity, 's' indicating the particle species, n the toroidal mode number, E s the energy of the single particle, n s the density, p s the pressure, e s the electric charge, r the minor radius coordinate, B the magnetic field, q(r) the safety factor, and the subscript 'h' refers to energetic ('hot') particles. It is worth noting that the ratio * w w h dh does not depend on the sign of the electric charge e s : thus, deeply trapped energetic electrons with a density profile peaked on-axis and of energy similar to that of energetic ions could be expected to drive a similar fishbone mode, propagating poloidally in the direction of the electron diamagnetic drift, i.e., opposite to the ion fishbone (although with some more unfavorable conditions [7], see section 3).
The first observation of fishbone oscillations driven by energetic electrons (e-fishbones) is reported almost two decades later in DIII-D [6]. In that experiment, strong MHD activity was observed in presence of NB ion heating, in conjunction with off-axis EC current drive and heating on high field side (HFS) and negative central shear equilibria with » q 1 min . The fishbone oscillations were stronger when EC was applied on the HFS equatorial plane (q p » res , with q res the resonant poloidal angle of the EC wave absorption location), and decreased while decreasing q res toward q p = 2 res . From the DIII-D experiment the following conclusions were derived: (1) from ray-tracing and Fokker-Planck calculations, it was shown that energetic electrons with hollow radial density profile were generated slightly internal to the = q 1 min surface, with a substantial fraction of barely trapped particles (i.e., particles, which spend most of their time near the tips of their banana orbits, and are preferentially heated when q p » res ); (2) the diamagnetic drift velocity of the energetic electrons (whose sign depends on  ( ) e p sign s s , with ( ) e sign s and p s the sign of the electric charge and pressure of species 's', respectively) is parallel to that of the on-axis peaked energetic ions produced by NBs; (3) the orbit averaged toroidal precession velocity (depending on ( ) e E sign s s ) of trapped energetic electrons, which is opposite to the one of the energetic ions for deeply trapped particles, reverses its sign when considering barely trapped particles [17,18], thus becoming parallel to that of deeply trapped energetic ions. As a conclusion, barely trapped energetic electrons with inverted radial density profile could meet the instability condition and drive a fishbone instability, in analogy with deeply trapped energetic ions with on-axis peaked radial density profile (here, the subscript 'Ee' stands for 'energetic electron'). Fishbone like fluctuations at higher frequency were also observed in COMPASS-D [19] driven by ECRH and LH.
Electron fishbones driven only by HFS off-axis (near q = 1) ECRH were observed also in HL-1M tokamak [20]; applying LH waves was found to enhance the fishbone, but only in conjunction with ECRH. Meanwhile, the observation of fishbone oscillations driven by only LH waves has been reported on FTU [7,21,22] and Tore Supra [23,24].
A careful experimental analysis to characterize the direction of propagation of the electron fishbone mode has been performed on HL-2A tokamak [25][26][27][28][29], where mode frequency of energetic electron driven fishbones during off-axis ECRH on both HFS and low field side (LFS) was studied. In particular, in [25] barely circulating energetic electrons were considered responsible to drive electron fishbones during off-axis ECRH on the LFS, whereas both barely circulating and barely trapped energetic electrons were considered responsible to drive electron fishbones during HFS heating. More recently, an interesting link between nonlinear dynamics of fishbone fluctuations in NB heated discharges and non-local thermal electron heat transport has been demonstrated in HL-2A [30].

Linear and nonlinear theory of e-fishbones
The theoretical framework for analyzing linear and nonlinear fishbone dynamics has been recently reviewed in [31]. Here, we briefly summarize the analysis given therein, referring to original works for more in-depth discussions.
Fishbone linear stability and nonlinear evolution can be described by the general fishbone like dispersion relation (GFLDR) [32,33] Here, we have assumed fishbone fluctuations with toroidal mode number n=1 and s denotes magnetic shear at the rational surface = r r s , where the safety factor = ( ) q r 1 s . For convenience, we adopt straight magnetic field line toroidal flux coordinates q z ( ) r, , , with r the radial (magnetic flux) coordinate, while θ and ζ are periodic angular variables in poloidal and toroidal directions, respectively. In equation (1), Λ accounts for the kinetic/ singular layer response at = r r s , where the fishbone mode structure is sharply varying. When s at = r r s vanishes, the singular layer response must be suitably rewritten into a form similar to equation (1), which can be found in [31][32][33] and is assumed to be adopted if necessary. Such explicit form is not given explicitly here, since it is not needed in the discussion of the formal properties of the GFLDR. Meanwhile, terms on the right-hand side of equation (1) describe potential energy contributions from the regular region, separating fluid (dŴ f ) and kinetic (dŴ k ) responses consistently with the original approach in [5]. Explicit expressions for Λ, dŴ f and dŴ k are given in [31][32][33] and will be omitted here for brevity.
The kinetic/singular layer response Λ accounts for the structures of the shear Alfvén wave continuous spectrum including kinetic and geometry effects, e.g., continuum damping, neoclassical inertia enhancement and Landau damping [32]. Meanwhile, dŴ f represents the potential energy fluctuation due to the fluid thermal plasma response and the fluid/convective behavior of the EP component. Kinetic responses of thermal and supra-thermal plasma are accounted for by dŴ k , which, e.g., describes instability drive by resonant wave particle interaction at . is taken along a closed equilibrium particle orbit.
Furthermore q q º ∮ ∮ q qd d [34] and ℓ indicates the 'bounce' harmonics of the considered wave-particle resonance.
In linear theory, equation (1) predicts that fishbones can be excited by magnetically trapped EPs at the precession frequency [5]. Thus, for deeply trapped EPs with pressure profile peaked at the magnetic axis, fishbones above excitation threshold are expected to rotate in the EP diamagnetic direction, as observed originally [4]. For e-fishbones, this situation is experimentally more difficult to achieve as they require a particularly strong fast-electron source, due to the stronger continuum damping for modes rotating in the electron than in the ion diamagnetic direction [7]. In fact, most favorable excitation conditions for e-fishbones are off-axis peaked EP pressure profiles with modes rotating in the ion diamagnetic direction, consistent with equations (2) and (3) and = ℓ 0 for barely trapped/barely circulating EPs affected by precession reversal [6]. The important role of finite -(¯) q 1 for barely circulating EPs in equation (3) has been recently emphasized in [35]. In the following sections, we will numerically investigate both situations; i.e., the case of pressure profile peaked on axis and e-fishbone excited by deeply trapped EPs, as well as the off-axis peaked EP pressure profile case excited by barely trapped/circulating electrons.
When nonlinear physics is investigated, equation (1) can be cast as a nonlinear equation for the evolution of the fishbone amplitude and essentially fixed (linear) radial mode structure [31]. The nonlinear wave-wave and wave-EP interactions are dominated, respectively, by the nonlinear responses L NL and dŴ k NL [31][32][33]. Both zonal flows and currents, i.e., in general zonal structures (ZSs), contribute to L NL and are generated by the dominant = ( ) ( ) m n , 1,1 component of the fishbone mode [36]. Meanwhile, dŴ k NL is due to phase space ZS (PSZS), that is self-consistent nonlinear distortions in the EP distribution function having the same symmetry (spatial structure) of the underlying EP reference equilibrium [34]. In general, dŴ k NL and nonlinear dynamics of PSZS dominate over L NL and ZS nonlinear response for sufficiently strong EP drive. However, the role of ZS becomes increasingly more important approaching the excitation threshold, and L NL has to be considered on the same footing as dŴ k NL near marginal stability [31]. This transition in the nonlinear behavior is quantitatively affected by kinetic and geometry effects, such as neoclassical inertia enhancement, which must be taken into account for a proper description of the self-consistent nonlinear fishbone dynamics in toroidal fusion plasmas. In this work, for simplicity, we consider a sufficiently strong supra-thermal electron source and numerically investigate only the effect of dŴ k NL and nonlinear dynamics of PSZS [34].
Following [31,34], and assuming deeply trapped resonant EPs with negligible orbit width and a rigid plasma displacement x d r0 , the evolution equation for the EP PSZS can be cast as (the magnetic moment) and  = v 2 2 (energy per unit mass) as well as r but, for conciseness of notation, only time dependence is indicated explicitly. Meanwhile, 1 0 i 0 denotes the Laplace transform. Furthermore, B 0 is the magnetic field on the magnetic axis at toroidal major radius = R R 0 and fishbone fluctuations are assumed to occur with the time dependent frequency w t g t + ( ) ( ) i 0 . The τ notation instead of t explicitly refers to sufficiently slow time variation, such that w gw  |˙| | | 0 0 [34]. Equation (4) coupled with equation (1) via dŴ k NL give the description of the self-consistent nonlinear evolution of PSZS and fishbone oscillations [31]. This process is generally non-perturbative; that is, fishbone spatiotemporal structures affect EP transport and vice versa. It is worthwhile noting that fishbone spatiotemporal structures are evolving in time even though the radial mode structure remains close to a rigid plasma displacement. In fact, mode amplitude and frequency are changing in time, and influence the phase space structure of ( ) F t 0 due to the radial and velocity space dependence of the resonance conditions, equations (2) and(3). With non-perturbative EP response, nonlinear evolution is dictated by maximization of wave-particle power transfer, as discussed in [31,34]. In particular, equation (4) suggests that PSZS evolution tries to preserve the resonance condition; i.e., that g t Q~- | | 2 NL 2 , with Θ the wave-particle phase and t NL the characteristic nonlinear time [31,34]. The ensuing frequency chirping is typically non-adiabatic; that is w g t~-|˙| 0 2 NL 2 and, thus, resonant EP continuously enter and leave the resonance region since wave-particle trapping is suppressed, consistently with maximization of wave-particle power transfer. A given group of particles remains in resonance for a finite interaction time, t I , given by the condition 2 ), resonant interaction is preserved during effective mode amplification and is rapidly lost after residual resonance detuning has shifted the wave-particle phase by p , similar to the slippage of resonant electrons with respect to a short free electron laser pulse [31,34]. In this way, mode amplification can continue as resonant EPs continuously enter (and leave) the resonance region. Equation (5) corresponds to a finite interaction length, Dr I , set by the distance traveled by the PSZS group velocity ( x w d r 0 0 ) over t I [31,34].
Strong fishbone excitation occurs when Dr r; I s i.e., when resonant particles are convectively pumped out from within the = ( ) q r 1 s magnetic flux surface as discussed in the original work by White et al [37]. In such a case, saturation corresponds to balancing the convective power loss with the wave-particle power transfer, . This behavior is demonstrated by present numerical simulation results with on axis peaked supra-thermal electron pressure profile, resonantly driven by deeply trapped particle at the precession resonance. For weaker EP drive and D < r r I s , after a group of particles looses the resonance condition according to equation (5), a new group of particles can enter the resonance region because of nonadiabatic frequency chirping ( w g t~-|˙| 0 2 NL 2 ) and mode amplification can continue until it is suppressed by equilibrium non-uniformity [31,34]. This behavior has been recently observed by numerical simulations based on a reduced description of the fishbone burst cycle that can be obtained from equations (1) and (4) with some additional simplifying assumptions [38].
When EP drive is further reduced, the effect of ZS should be considered on the same footing of PSZS nonlinear dynamics, as anticipated above. This is beyond the scope of the present work, where we focus on nonlinear wave-EP interaction and the effect of dŴ k NL only. Application of equations (1) and (4) in such a case, approaching marginal stability and with perturbative EP drive, yields the prediction fishbone saturation by local EP redistribution. In the case of supra-thermal electron pressure profile peaked off-axis, numerically analyzed in this work as illustration of e-fishbones driven by barely trapped/circulating EPs, w res radial profile has a local maximum at r s , corresponding to slowly evolving (adiabatic) PSZS of 1 3 radial extension, with w res0 the second radial derivative of w res at r s . It is then possible to show that fishbone saturation due to local EP relaxation occurs at

Description of XHMGC
In the following sections we will present the results of numerical simulations performed using the HMGC code [39][40][41], which is a hybrid [42] MHD gyrokinetic code originally developed at the Frascati Laboratories. In HMGC, the thermal plasma is described by  ( ) O 3 nonlinear reduced MHD equations [43], which describe circular shifted magnetic surface equilibria; moreover, the limit of zero bulk plasma pressure is also assumed; the energetic particles are described by nonlinear Vlasov equation in the drift-kinetic limit, solved using particle-incell technique, the two components (thermal and energetic particles) being coupled [42] via the pressure tensor term of the EP species entering in the extended momentum equation of the bulk plasma. The hybrid scheme allows to consider the effect of the energetic particles on the electromagnetic fields self-consistently, i.e., they are retained non-perturbatively: thus, the mutual effect of energetic particles and MHD-like modes (as, e.g., toroidal Alfvén modes, TAEs and internal kink modes), as well as modes which do not have their MHD counterpart (as, e.g., energetic particle driven modes, EPMs), can be studied properly; in particular, energetic particles will contribute to both time evolution of the mode (i.e., to the growth rate and frequency) and to its spatial structure (i.e., to the eigenfunction). HMGC has been extensively used to study TAEs and EPMs [41,44,45], as well as in the analysis of modes observed in existing devices (JT-60U [46], DIII-D [47]) or expected in forthcoming burning plasmas (ITER [48,49]) and proposed experiments (FAST [50,51]). The original version of HMGC has been recently extended to include new physics (XHMGC [9]): diamagnetic effects and thermal ion compressibility are retained in the extended momentum equation of the bulk plasma through the divergence of the thermal ion pressure tensor, obtained by solving the nonlinear Vlasov equation for that population, in order to account for enhanced inertia response [7,10,11] and ion Landau damping [12]. Moreover, XHMGC is able to treat simultaneously, using the kinetic formalism, up to three independent particle populations, assuming different equilibrium distribution functions (as, e.g., bulk ions and electrons, energetic ions and/or electrons accelerated by NB, ICRH, ECRH, fusion alphas, etc). The XHMGC code has been also used to simulate fishbone modes driven by energetic electrons [8]. With respect to energetic ion driven modes (as, e.g., TAEs and EPMs), the simulation of e-fishbones poses the challenge of properly follow the extremely fast parallel electron motion along equilibrium magnetic field lines, in order to correctly evaluate the effective bounce averaging during both linear and nonlinear dynamics. To this purpose, a suited sub-cycling scheme has been introduced and solved in the particle-in-cell scheme. As synthetic diagnostic tool, XHMGC allows to follow, in a self-consistent simulation, a set of test particles; the phase-space coordinates of such particles are stored in time, and can be used to compute a variety of single particle physical quantities as , e.g., the single particle frequencies of the suprathermal electrons, namely, the precession and bounce frequencies. The resonances underlying the linear instability can be clearly identified in this way. Furthermore, the use of EP phase-space diagnostics, based on Hamiltonian mapping techniques [13] generating kinetic Poincaré plots, allow us to isolate the physics processes underlying fishbone mode saturation, frequency chirping and secular (versus diffusive) EP redistribution.
The energetic electrons ('Ee') distribution function used in the following simulations is: In the following simulations performed with XHMGC, the contribution to Landau damping, enhanced plasma inertia (mostly due to trapped particles), finite compressibility of thermal ions, as well as Landau damping and finite compressibility of thermal electrons will be all treated kinetically by considering isotropic being the corresponding density and temperature profiles, with = j i e , . We will neglect mode-mode coupling nonlinearities, thus considering single n toroidal mode number simulations, while particles nonlinearities will be fully retained.

Energetic electrons with density profile peaked on-axis
As a first example of e-fishbone we will consider an energetic electron population with on axis peaked density profile. Similarly to the conventional energetic ion driven fishbones, deeply trapped energetic particles are expected to drive the mode. The same FTU-like equilibrium and scenario of [8] will be considered in this section, namely a torus of circular cross section with inverse aspect ratio  = » a R 0.3 0 (with a and R 0 the minor and major radius, respectively). The safety factor profile is slightly reversed, with »  , while the onaxis bulk electron Larmor radius is r »´a 1.32 10 e 4 . Energetic electrons described by a distribution function with perpendicular temperature much higher than the parallel one are considered in the simulations by assuming a = cos 0 0 and D = 0.1, see equation (7), in order to maximize the fraction of trapped particles, which are expected to contribute to the resonant excitation of electron fishbone mode; the energetic electron density radial profile is , whereas the energetic electron temperature is assumed to be uniform , while the on-axis energetic electrons Larmor radius is r »´a 3.52 10 Ee 4 . This equilibrium has been already analyzed in [8], were it was shown that the electron fishbone mode was destabilized above a certain threshold energetic electron density, propagating poloidally in the direction of the energetic electron diamagnetic velocity (which is, for this equilibrium, also parallel to the bulk electron diamagnetic velocity), and excited by resonance with deeply trapped energetic electrons. Here, we will reconsider the linear results presented in [8], where the kinetic contribution of the energetic electrons and bulk ion was considered, by also adding the kinetic contribution of bulk electrons. Moreover, a novel nonlinear analysis using test particles Hamiltonian mapping (TPHM) techniques [13,14] will allow us to illuminate the nonlinear saturation and radial transport associated with such mode.

Linear dynamics
In this section we will investigate the relative importance of different driving and damping processes accounted for in the model, i.e., energetic electrons compressibility, thermal ion compressibility and diamagnetic effects, and thermal electron compressibility. Following [9], where the model implemented in XHMGC has been described in detail, let us consider the perpendicular component of the extended MHD momentum equation:  figure 2, corresponding to switching on, one after the other, the contributions isolated in equation (8). First, only the divergence of the energetic electron pressure tensor ( · ) P Ee is retained, then also the diamagnetic bulk ion term is added and, subsequently, the divergence of the thermal ion pressure tensor ( · ) P i , which account for the thermal ion Landau damping and generalized inertia, retaining consistently the actual dynamic response of trapped and circulating thermal ions (see also section 2.2 and appendix A of [7]). Finally, the divergence of the thermal electron pressure tensor ( · ) P e is also included. The contribution of energetic electrons drives the mode, which has a clear internal kink characteristic with a dominant m=1 component localized, in radius, approximately inside the q min surface ion term, very little variation is observed, both in growth rate and frequency: indeed, the absolute value of this term, evaluated at its maximum radial position ( » r a 0.35) is much less (by a factor »30) than the absolute value of the frequency of the mode. When adding the term ( · ) P i , on the contrary, the growth rate of the mode is notably reduced, showing as the effect of considering the thermal ion Landau damping and enhanced inertia increases the threshold in n n Ee0 i0 required to destabilize the mode; also the absolute value of the frequency of the mode increases. Finally, when adding the term ( · ) P e which accounts for the bulk electrons, an increase of the growth rate is observed, which diminishes its importance as n n Ee0 i0 is increased. Note that the positive contribution to the growth rate of the bulk electrons appears only if the electron fishbone is already driven unstable, indeed in absence of the energetic electrons the system is still stable.
Some more insight on the linear dynamics of the electron fishbone driven by energetic electrons with density profile peaked on-axis can be gained considering the power exchange between the various particle species and the wave. In figure 4 the power exchange in the (v M , ) space is shown, for several radial shells, in a simulation in which the kinetic contributions of all the particle species are retained. The dominant drive contribution comes from the deeply trapped energetic electrons (the solid and dashed curves in each plot refer to the approximative boundary between trapped and circulating region in the phase space, for the inner (solid curve) and outer (dashed curve) radii of the radial shell considered, respectively). Note that in the plots shown in figure 4 the normalized parallel velocity u refers to the parallel velocity the simulation markers ('macro-particles') have when crossing the equatorial plane at the poloidal angle q = 0; in particular, for trapped particles, the external 'leg' of the 'banana' orbit is chosen. Let us first consider the power exchange relative to the energetic electrons: as already stated in [8], the drive is given by the deeply trapped electrons, figure 4, upper-left plot. Making use of the test particle Hamiltonian mapping techniques [14], we can gain some insight on the characteristic resonances of the energetic electrons. As stated in section 4.1, the code HMGC can be used to evolve a set of test particles in the  time dependent e.m. fields computed by a self-consistent simulation [14]. A suitable choice will be to choose the set of test particles as the ones with strongest resonance with the wave during linear growth phase. This can be done by choosing the particles which belong to the phase space region where the maximum power exchange between energetic particles and wave occurs (see, e.g., figure 4). Once the coordinates of such a region have been selected (i.e., = r r 0 , = M M 0 , =  v v ,0 ), we can define, following [13,14], the quantity w º - being the canonical toroidal angular momentum. The quantity C is a constant of the (perturbed) motion, provided that the perturbed field is characterized by a single toroidal mode number n and a constant frequency. At the leading order [14], we can approximate , for n = 1, = ℓ 0 (with ℓdenoting the 'bounce harmonics'), is satisfied in correspondence of the peak of the power exchange between the (test) energetic electrons and wave, confirming that the precession frequency w d of deeply trapped energetic electrons is driving the mode. Note also that the radial profile of the test particles resonant frequency w ( ) r res , after a strong variation with radius for  r a 0.15, is quite flat and close the value of the (linear phase) mode frequency w 0 : this feature will be relevant for the evolution of the mode during nonlinear saturation (see section 4.2.2).
From figure 4 the bulk ions are shown to contribute to the damping of the mode, mainly with co-and counter-passing particles, whereas the bulk electrons, depending on the radial shell considered, contribute to damping the mode with both passing and trapped particles in the internal radial shells, but strengthen the drive some where outside the q min radius. Indeed, applying the same TPHM technique used above, it is possible to show that, in correspondence with the positive contribution of the bulk electrons in the radial shell » r a 0.46, which exhibit a maximum for the power exchange at

Nonlinear dynamics
In the present section we will consider the nonlinear dynamics of the electron fishbone driven by energetic electrons with density profile peaked on-axis. We will refer, in this section, to a set of simulations in which only the contribution of the thermal ions treated kinetically will be considered beside the one of the energetic electrons, neglecting thus the kinetic response by thermal electrons, as in the case of the simulations presented in   (the weakest simulation shown in figure 2).
As already discussed in [8], the saturation is characterized by a pronounced downward (in absolute value) frequency chirping, as it is shown in figure 8, where the frequency spectra of the electrostatic potential j w ( ) r, are shown at several times of the simulation (during linear phase ( w » t 400

A0
), end of the earlier phase of )). In figure 8, superimposed on the frequency spectra, the instantaneous resonant frequency w w = ( ) r t C M , ; , res res 0 0 is also reported: it is worthwhile noting that the instantaneous resonant frequency w ( ) r t , res changes in time, being always close to the local lower Alfvén continuum, and that the radial extension of the power spectrum j w | ( )| r, 2 is almost unchanged during the downward (in absolute value) chirping, as a consequence of the stiffness of radial shape of the internal kink eigenfunction. Furthermore, the mode, although chirping down, is always superimposed to the instantaneous w ( ) r t , res curve, thus suggesting phase-locking with the energetic particles identified by = = C C M M , 0 0 . In figure 9 the spectrogram of the mode is shown at the radial location where the electrostatic potential has its peak: the frequency chirps down (in absolute value) as the simulation leaves the exponential growth of the linear phase (for  w t 600

A0
), in correspondence to bursts observed in the total MHD energy. By following the same set of test particles considered in section 4.2.1 also in the nonlinear phase of the simulation, we can further investigate the saturation mechanism of the mode and the associated radial transport of resonant particles. In figure 10 (left) the radial density profiles of the test particles are shown, for several times during the simulation ( w = t 400

A0
, during the saturation phase). Radial flattening is strongly evident in these resonant particles, starting at the radial position where the power exchange between energetic particles and wave is stronger ( » r a 0.185, see figure 5). It is worthwhile noting that this strong radial transport of the energetic particles can become much less evident after averaging (summing) over the total (resonant plus non-resonant) EP radial density profile. From figure 10 (left) it can also be observed that the flattening of the radial profile, that starts at time w » t 600 , as computed from the test particles evolution, continues to flatten in radius with respect to the linear phase and extends outward, but also, because of the frequency chirping of the mode, it shifts up: as a consequence, the resonant radius moves outward, i.e., toward higher values of the canonical toroidal angular momentum f P . In figure 11 the kinetic Poincaré plots are shown for the same times of simulation of figure 10 with test particles colored accordingly to their initial f P value: here, whenever a test particle completes a full banana or transit orbit in the poloidal plane at =t t , by crossing the outer equatorial plane (q = 0), the corresponding values of the canonical toroidal angular momentum f (ˆ) P t and the wave-particle phase Q(ˆ) t are computed, with and f being the toroidal angle. Each test particle position in the plane (Q f P , ) is then updated by moving the marker in the 'kinetic Poincaré' plot [13]. The extremes of the vertical axis = - 0.13, 50 correspond to,  , the resonant particles that are displaced at larger f P (because ofÉ B drift), will, in turn, continue to drift outward (instead of drifting horizontally in Θ, and reversing their motion toward the resonant layer when feeling a change of sign of theÉ B drift). This outward drift gives rise (see figure 11, w = t 700, 800, 900, 996

A0
) to vertical elongated structures, i.e., to a large radial transport of energetic resonant particles (see also the almost flat test particle density profiles for the same times in figure 10 (left)), consistent with the theoretical analysis given in section 3 [31,34]. This process ends when the flattening of the test particle density profile approaches the q min radius, where the internal kink eigenfunction sharply decreases ('radial decoupling', see [14]); it has to be noted that very little variation is observed in time on the shape of the eigenfunction, being the internal kink type eigenfunction quite stiff.
Evidence of phase locking is shown in figure 12, where the average of the precession frequency w ( ) t d of the linearly resonant particles, weighted with their power transfer during linear phase is shown, together with the width of its distribution (see [8]), and compared with the time varying (chirping) frequency of the mode (see figure 9). Also, the mode adjusts its frequency (thus keeping Q » |˙| 0, the 'phase-locking' condition) in order to remain 'tuned' with the resonant particles which, meanwhile, experience their outward displacement.
Finally, we plot in figure 13 the saturation amplitude of the = = m n 1, 1 Fourier component of the electrostatic potential j sat 1,1 as the strength of the EP linear drive varies. For the saturation amplitude, we took Figure 11. Kinetic Poincaré plots with test particles colored corresponding to their birth value of f P . Note that the extension of the Θ axis has been doubled and test particle markers has been replicated in the domain  p p Q < 2 4 to enhance the readability of the plots. The arrows in the first plot indicate the direction of the particle drift along Θ above and below the resonant layer º » the maximum value of j 1,1 at the first plateau in time (see, e.g., figure 7, where the first plateau occurs at w » t 720

A0
). From the simulations we can infer that j g w µ a | | ( | |) . These results compare favorably, for weak drive, with the findings of [52], whereas, for sufficiently strong drive, are in fair agreement with that given in [31] and already anticipated in section 3

Energetic electrons with density profile peaked off-axis
In this section the first global hybrid MHD-Gyrokinetic simulations of electron fishbones driven by energetic particles with density profile peaked off-axis [53] will be presented. This kind of equilibria is closely related to the experimental configuration in which electron fishbones have been observed in current devices. In these experiments, HFS off-axis heating is applied close to the q min flux surface in the equatorial plane, using ECRH; thus, an inverted (positive) gradient of the energetic electron density profile is generated in the radial region of the discharge which is internal to the q min flux surface and in which the internal kink can develop. Moreover, because of the HFS deposition, a selective heating on barely trapped/circulating particles will be obtained [6].
Recalling the stability condition, * w w > 0 Ee [5], and noting that * w Ee depends on the sign of the radial gradient of the energetic electron pressure profile, instability can occur only by resonance with energetic electrons characterized by precession reversal; i.e., barely trapped/circulating energetic particles [17,18]. The equilibrium considered here has the same bulk density and temperature profiles and plasma parameters of the peaked on-axis one (see section 4.2 and figure 14 right), except for the inverse aspect ratio,  = 0.1, and the  safety factor profile q, which also in this case is slightly reversed with » q 1. has been used in order to minimize the continuum damping and facilitate the occurrence of the energetic electron driven fishbone [7,35]. Moreover, safety factor profiles with a reversed shear is known to enhance the reversal of precessional drift [17,18,54].
The energetic electrons radial density profile, peaked off-axis, has been chosen as y = as for the energetic electron density profile peaked on-axis, see equation (7), and thus the distribution function is still symmetric in velocity space, and no net current is driven by energetic particles. The choice of D = 0.5 is such to ensure the presence in the energetic electrons distribution function of a sufficient fraction of barely trapped/circulating particles (see figure 15). In the following simulations the toroidal mode number is n=1, and the poloidal Fourier components retained are = ¼ m 1, , 6, normalized resistivity =-S 3 10 1 5 and viscosity nt =´a 3 10 A0 2 8 , as before. Besides the kinetic contribution of the energetic electrons ( · ) P Ee , both the diamagnetic bulk ion term, and the thermal ion compressibility ( · ) P i treated kinetically are retained, whereas the contribution ( · ) P e of the bulk electrons will be neglected for simplicity, see equation (8). Moreover, the choice of a shaped energetic electron  temperature profile, which strongly decreases outside the » q 1 surface, has the beneficial effect of inhibiting the growth of modes with dominant poloidal mode numbers higher than unity, which can be driven unstable by deeply trapped energetic electrons outside the q=1 surface, where the energetic electron density gradient becomes negative.

Linear dynamics
In this section we present linear dynamics results for the electron fishbone driven by energetic electron with offaxis peaked density profile. In figure 16 the radial profiles of the Fourier components, the poloidal structure and the power spectrum of the electrostatic potential are shown, for = n n 0.0095 Ee0 i0 : the radial structure of the poloidal Fourier components is dominated by the m, n=1,1 component, which is localized in the region  q q min , showing the characteristic shape of the internal kink radial displacement (x j µ r r ). The structure of the mode in the poloidal plane rotates in time in the clock wise direction, i.e., opposite to the direction observed in the peaked on-axis EP radial profile (see section 4.2), which corresponds to a mode propagating in the direction of the diamagnetic velocity of the energetic electrons (which is parallel, for a peaked off-axis energetic electrons density profile, to the direction of the bulk ion diamagnetic velocity). The frequency of the mode is quite low, see figure 16 (right), as expected for equilibria with low values of Dq (see, e.g., [35]). In figure 17 the growth rate and frequency of the electron fishbone driven by energetic electrons with radial density profile peaked off-axis is shown: a linear dependence on n Ee0 is observed, above the threshold » n n 0.007 Ee0 i0 , with a real frequency weakly dependent on n Ee0 .
From the power exchange between the energetic electrons and wave, as shown in figure 18, it is possible to infer which fraction of energetic particles is driving the mode. In figure 18 (left) the power exchange in the radial shell   r a 0.33 0.41 is shown, with the curves approximating the trapped/untrapped boundary (solid black) and the barely/well circulating boundary (solid red) superimposed. Here, we follow the definition given in [7] for the barely circulating particles, defined as , with:   Indeed, the maximum power exchange occurs for particles in the region of velocity space belonging to that of barely circulating ones (in particular the counter-circulating ones, red pattern), with some minor contributions coming from the well circulating particles, both co-and counter-circulating (outside the solid red curve, green pattern); trapped particles, on the contrary, give a damping contribution (light blue to dark blue patterns) to the mode, as expected (see figure 18 (right), where only the power exchange due to trapped particles (k > 1 2 ) is shown in order to enhance the relative size of their contribution).
An analysis, similar to the one shown for the peaked on-axis energetic electrons density profile case, has been done for this simulation, in order to identify the characteristic resonance responsible for driving the mode. In particular, the Hamiltonian mapping technique has been applied to a set of test particles defined by the C M ,  . Contribution from the full population in the radial shell   r a 0.33 0.41 (left), and only trapped particles (right), same radial shell; black lines refer to the boundary between trapped/untrapped particles, whereas solid red curve refers to the boundary between barely circulating and well circulating ones (for  r a 0.41). . Over-imposed is the radial dependence of the test particle power exchange, (blue dotted-dashed curve). radial locations ('double resonance'), as a consequence of the -( ) q 1 term in the resonant condition for circulating particles. Note that the variation of w ( ) r res around the resonant positions is stronger than the one observed in the peaked on-axis EP density profile (being, for the two cases considered, res,off axis res,on axis ); the effects of this difference will be manifest during the nonlinear saturation phase w.r.t. the peaked on-axis case.

Nonlinear dynamics
The evolution of the simulation with = n n 0.0095

Ee0 i0
is shown in figure 20, where the time evolution of the total (kinetic plus magnetic) volume integrated MHD energies, W tot;m,n , for the m n , Fourier components considered in the simulation are shown.
In figure 21  ). Note that the double resonance structure, which has the sign of w ¶ ¶r res at one resonant radius the opposite of that at the other resonant radius, makes particles in resonance at the inner radius rotating in the opposite direction w.r.t. the ones in resonance at the outer radius (see figure 21, second frame from left). While entering the fully nonlinear phase ( , see figure 21, third frame from left), we note that the two island structures tend to insinuate oneself into the other, having a f P extension (or, equivalently, a radial extension) of the order of the distance between the two resonance layers - ,res1 . As the test particles are displaced outside the resonant layer (toward < r r res1 , or > r r res2 ), where the characteristic resonant frequency changes rapidly with radius (see figure 19), even changing its sign and, thus, not satisfying any more the instability condition * w w > 0 Ee , the mode has no 'convenience' in adjusting its frequency to that of the linearly resonant particles. Indeed, little variation of the frequency during the saturation phase is observed, , in the plane Q f ( ) P , . Test particles are colored according to their initial f P values. The arrows in the first plot indicate the direction of the particle drift along Θ above, in between, and below the resonant layers » f P 7 ,res . and the saturation of this simulation can be ascribed to 'resonance detuning' (see, e.g., [14,31,34,55]). The result of saturation on the radial density profile of the energetic particles is shown in figure 22: a flattening of the profile in a region somewhat larger than the one between the two resonant surfaces » r a 0 , see also the Poincaré plots, figure 21, third and fourth frame from left).
In figure 23 the scaling of the saturation amplitude of the electrostatic potential versus the ratio of the linear growth rate to the frequency of the mode g w L 0 is shown, in a scan in which the EP density is varied: a stronger scaling, j g w » ( )

Conclusions
In present days, e-fishbones have been observed on a variety of devices: they are internal kink-type instabilities driven by energetic electrons generated by, e.g., EC and/or LHH/current drive. A brief summary of experimental evidences have been given in section 2. Similarly to the well known ion fishbones, e-fishbones are  driven by wave-particle interaction at magnetic toroidal precession frequency of energetic electrons: they are characterized by linear dispersion relation similar to that of ion-fishbones (see section 3 for a brief summary of the theory of linear and nonlinear dynamics of e-fishbones). The large interest raised about the e-fishbones is also related to the fact that the radial transport, induced by fluctuations driven by magnetically trapped resonant particles, due to precession resonance, depends on energy and not mass of the particles themselves. Also, the energetic electrons on present devices are characterized by having much smaller resonant particle orbit width w. r.t. the characteristic scale length of the fishbone instability (which is of the order of the rigid radial displacement of the internal kink-like eigenfunction) similarly to what is expected for energetic ions (e.g., fusion α-particles) in burning plasmas, because of large plasma current, and differently from what can be obtained for energetic ions in present experiments. Thus, e-fishbones offer the opportunity to study both experimentally and numerically linear and nonlinear dynamics of wave-particles interactions relevant for fusion plasmas also on present, not ignited devices, and to compare with theory.
In the present paper, linear and nonlinear numerical simulations of fishbones driven by energetic electrons have been presented, using the hybrid MHD-Gyrokinetic code XHMGC [9]. Equilibria with both on-axis [8] and off-axis [53] peaked energetic electron radial density profile have been considered. For the on-axis peaked energetic electron radial density profile case, it has been enlightened the effects of considering the kinetic thermal ion compressibility and diamagnetic effects (in order to allow for an entirely novel treatment of enhanced inertia response and ion Landau damping) and the kinetic thermal electrons compressibility effects. The most important effect comes from the thermal ion compressibility treated kinetically, resulting in a considerable enhancing of the instability threshold energetic electron density. Resonances between wave and energetic electrons and thermal species have been analyzed using TPHM techniques [13,14], clearly identifying the precession frequency w d of deeply trapped energetic electrons in driving the mode [8]. Nonlinearly, the mode experiences a strong downward frequency chirping, adjusting its frequency in order to remain tuned with the resonant particles which are displaced outwardly; this results in keeping Q » |˙| 0, the phase-locking condition, and induces the simultaneous flattening of the radial density profile of the resonant energetic particles together with the outward displacement of the resonance radial location. The process ends when the flattening of the resonant particles density profile approaches the q min radius, where the internal kink eigenfunction, that remains almost unchanged in its shape during nonlinear phase, sharply decreases in amplitude ('radial decoupling', see [14]). Note that some of these features (chirping frequency, ejection of resonant particles) are consistent with previous findings of ion fishbone simulations [56]. The nonlinear saturation of the dominant Fourier component ( = m n 1 1) of the electrostatic potential, j sat1,1 , has been shown to scale almost proportional to the strength of the ratio of the linear growth rate to the mode frequency, j g w µ sat1,1 L 0 , for sufficiently strong drive (  g w 0.15 L 0 , for this particular simulation case), whereas stronger dependence is obtained in the weaker drive regime. The scaling obtained for strong drive, when the finite interaction length of the phase space zonal structures is of the order of the q min radius, r s , is consistent with the analytical findings [31], as reported in section 3.
In this paper we have also presented the first global hybrid MHD-Gyrokinetic simulations of electron fishbones driven by energetic particles with density profile peaked off-axis. Equilibria with off-axis peaked energetic electrons radial density profile are more closely related to the experimental conditions in which e-fishbones have been observed in current devices, where HFS off-axis heating using, e.g., ECRH forms barely trapped and/or barely circulating energetic electrons. These particles, because of the the toroidal precession reversal of their motion, can fulfill the instability condition * w w > 0; Ee the mode is observed to rotate poloidally in the direction of the diamagnetic velocity of the energetic electrons (which is parallel, for a peaked off-axis energetic electrons density profile, to the direction of the bulk ion diamagnetic velocity). In this equilibrium the radial structure of the resonance for the energetic electrons driving the mode is dominated by the -( ) q 1 term entering in the expression of the resonance of circulating particles, see equation (3), thus having a double resonance in the vicinity of the q min radial location r s . Nonlinearly, kinetic Poincaré plots show a double island structure, one rotating in the opposite direction of the other, and shifted by π along the waveparticle phase Θ. Local flattening of the radial density profile of the resonant particles is observed in the radial positions of the double resonance, and no major frequency chirping is observed during saturation, which can be attributed to resonance detuning [14,55]. The scaling of the saturated amplitude, for this equilibrium, is in good agreement with the predictions given in section 3: j g w » ( ) sat1,1 L 0 3 for weak drive, and j g w » ( ) sat1,1 L 0 3 2 when w w w g D º -| | res0 0 L . It is worth noting that in this paper we have neglected the role of MHD nonlinearities, which are known to increase their importance near marginal stability [16,36]. Thus, a more complete simulation study will be required to explore in detail the nonlinear dynamics also in these regimes of weak energetic electron drive.