Toroidal Alfvén eigenmode triggered by trapped anisotropic energetic particles in a toroidal resistive plasma with free boundary

The toroidal Alfvén eigenmode (TAE), excited by trapped energetic particles (EPs), is numerically investigated in a tokamak plasma, using the non-perturbative magnetohydrodynamic-kinetic hybrid formulation based MARS-K code (Liu et al 2008 Phys. Plasmas 15 112503). Compared with the fixed boundary condition at the plasma edge, a free boundary enhances the critical value of the EPs kinetic contribution for driving the TAE. Free boundary also induces finite perturbations at the plasma edge as expected. An anisotropic distribution of EPs, in the particle pitch angle space, strongly enhances the instability and results in a more global mode structure, compared with the isotropic case. The plasma resistivity is also found to play a role in the EPs-destabilized TAE. In particular, the mode stability domain is mapped out, in the 2D parameter space of the plasma resistivity and a quantity defining the width of the particle distribution in pitch angle (for anisotropic distribution). A resonance layer in the poloidal mode structure, with the layer width increasing with the plasma resistivity, appears at the large width of the particle distribution in pitch angle space. A mode conversion, from the modified ideal kink by the EPs kinetic effect to the TAE, is also observed while increasing the birth energy of EPs. Computational results suggest that the TAE mode structure can be modified by certain key plasma parameters, such as the EPs kinetic contribution, the equilibrium pressure, the plasma resistivity, the distribution of EPs, as well as the birth energy of EPs. Such modification of the eigenmode structure can only be obtained following the non-perturbative hybrid approach (Wang et al 2013 Phys. Rev. Lett. 111 145003, Wang et al 2015 Phys. Plasmas 22 022509), as adopted in this study. More importantly, numerical results show that near the marginal stability point, the dominant poloidal harmonics of the TAE overlap with each other, and are localized at the tip positions of the Alfvén continua. This kind of TAE structure in high beta plasma with unstable ideal kink is substantially different from that of the conventional TAE.


Introduction
The toroidal Alfvén eigenmode (TAE) exists in the gaps of the Alfvén continuum spectra, which are formed by toroidal coupling of Alfvén waves with adjacent poloidal mode numbers in toroidal plasmas [1]. TAE is often easily excited by energetic particles (EPs), generated either by fusion products or by auxiliary heating such as the neutral beam injection (NBI) and the ion cyclotron resonant heating (ICRH) [2][3][4][5][6][7][8][9]. This kind of global instability can induce the redistribution and/or loss of EPs, with the latter having a direct consequence of increasing the heat flux on the plasma facing materials [10,11]. In fact, the TAE instability is considered to be one of the key factors degrading the alpha particles confinement in future devices such as ITER [12]. The transport and the loss of EPs, induced by TAE, have been computationally predicted [13][14][15][16][17][18][19] and experimentally observed [5,6,20]. Whilst the stability of the TAE is largely determined by the balance between the drive produced by EPs and continuum damping [21,22], and by other effects such as the resistive damping [23], the self-consistent modification of the mode structure, by various plasma parameters, has been less explored. Detailed understanding of the TAE structure, however, is of significant importance in terms of its nonlinear consequences, e.g. the EPs confinement.
Conventionally, the TAE study assumes a fixed boundary condition at the plasma surface [4,[24][25][26]. This boundary condition forces the radial plasma displacement to vanish at the plasma edge. This unrealistic assumption, neglecting the vacuum gap between the plasma and a conducting wall, may be acceptable for higher-n TAE (n is the toroidal mode number), but can be problematic for low-n modes, in particular for the n = 1 TAE, which is the focus of the present study. The n = 1 TAE instability is observed in both the conventional [27] and the spherical tokamaks [28].
In this work, we consider a free boundary condition, by which we mean the presence of a vacuum region between the plasma surface and a conducting wall, with the wall possessing a finite conductivity. The effect of free boundary condition on the TAE mode was numerically investigated by Chen et al [29], showing that the TAE frequency decreases with increasing the radial location of the conducting wall further away from the plasma. But their study did not consider the kinetic contribution from trapped EPs.
In this work, we shall apply the full toroidal stability code MARS-K to systematically investigate the n = 1 TAE excited by trapped energetic particles. Being a long-wavelength instability, the n = 1 TAE can produce significant energetic particle re-distribution, or even induce large scale energetic particle losses through profile avalanche. For instance, previous work showed that the dominant loss of counter-passing particles in ITER is produced by global low n TAE modes, predominantly the large n = 1 TAE [30]. On the other hand, TAE modes are mainly produced by trapped energetic particles. For instance, Zheng et al [9] showed that the precessional drift resonance is the dominant one for the TAE instability excited by trapped energetic ions. Both the free boundary condition and contributions from EPs are considered. The non-perturbative MHDkinetic hybrid approach is adopted [31], which allows a self-consistent modification of the eigenmode structure due to the presence of EPs [32]. Our computational results show that the free boundary condition has a stabilizing effect on the TAE, thus enhancing the critical value of the EPs kinetic effect for triggering the TAE instability, as compared to the case with a fixed boundary condition. We find that the TAE structure changes, as the mode eigenvalue varies as a function of the EPs kinetic contribution or the bulk plasma pressure. More specifically, the mode structure becomes more localized as the EPs kinetic effect decreases. The dominant poloidal harmonics overlap with each other near the marginal stability point for the mode, resulting in a mode localization near the tip positions of the Alfvén continua. The peak location of mode structure along the minor radius shifts inward with increasing the bulk plasma pressure.
Anisotropy of the EPs equilibrium distribution in the particle pitch angle space also affects the TAE growth rate. Compared with the isotropic case (such as that for α particles), the TAE is found to be more unstable as anisotropy increases. In contrast with previous study showing that the TAE is still unstable with the resistive plasma model [23], our systematic scans find a stability domain in the 2D parameter space of the plasma resistivity versus the width of the particle distribution in pitch angle. The critical value of the plasma resistivity, for the TAE stability, decreases as width of the particle distribution in pitch angle increases. The width of the resonance layer of the poloidal mode structure, induced by continuum damping at wide distribution, increases with the plasma resistivity. Finally, our numerical results show a smooth transition from an unstable ideal kink, which is modified by the EPs kinetic effect, to the TAE, with increasing EPs birth energy. Such a study mostly represents theoretical interests.
The paper is organized as follows. Section 2 describes the MHD-kinetic hybrid model implemented in MARS-K. The first part of section 3 presents the equilibrium. The second part of section 3 reports key novel numerical findings, including the effect of the boundary condition, the plasma pressure, the width of particle pitch angle for the anisotropic distribution, as well as the plasma resistivity, on the TAE instability. The third part of section 3 reports the other interesting results from systematic parameter scans, such as for the wall position and the birth energy of trapped energetic ions. Summary and a brief discussion follow in section 4.

Computational models
The MARS-K code numerically solves the linearized, single fluid MHD equations, with self-consistent inclusion of drift kinetic effects from multiple particle species [31]. This code has indeed been extensively benchmarked against other codes for the resistive wall mode instability, including MISK, PENT codes [33,34], and HAGIS [35]. All the benchmarking includes the EPs kinetic effect. The core equations of the MARS-K model are where the variables ξ, υ, Q, j, p represent the plasma displacement, the perturbed velocity, magnetic field, current and pressure tensor, respectively. ω is the complex eigen-frequency of the mode, corrected by a Doppler shifted frequency nΩ, with n being the toroidal mode number, Ω the plasma rotation frequency along the toroidal angle φ. ρ is the plasma equilibrium density, and η is the plasma resistivity. The equilibrium magnetic field and plasma current density are denoted by B and J, respectively. R is the plasma major radius. Z is the unit vector in the vertical direction, and I the unit tensor. b = B/B, B = |B| . The kinetic terms enter into the MHD equations via the perturbed parallel and perpendicular pressure components, p (ξ ⊥ ) and p ⊥ (ξ ⊥ ) which depend on the mode eigenfunction as well as eigenvalue, as shown below. The perturbed pressure tensor p is self-consistently included into the model via the momentum equation (2). With an assumption that the perturbation varies as exp(−iωt + inφ), the perturbed pressure can be written as here σ =⊥, , E = Mυ 2 /2, and E ⊥ = Mυ 2 ⊥ /2. M is the particle mass, υ and υ ⊥ are, respectively, the parallel and perpend icular (to the equilibrium magnetic field) velocities of the particle guiding center motion. The perturbed distribution function of particle species j is denoted by f 1 j , with j = i, e, h labeling the thermal ions, thermal electrons, and energetic particles, respectively. The integral is carried out over the particle velocity space Γ. Each component of the perturbed pressure involves both adiabatic (superscript 'a') and nonadiabatic (superscript 'na') parts: The perturbed distribution function, defined in the Lagrangian frame, satisfies the following equation where f 0 j (ε, ψ) is the equilibrium distribution function of particles, which we assume to be Maxwellian for thermal ions and electrons, and slowing down distribution for energetic ions. f 0 ε and f 0 P φ denote the derivatives with respect to the particles energy and toroidal canonical momentum, respectively. ψ is the poloidal flux, ε = ε k + ZeΦ = Mυ 2 /2 + ZeΦ is the particle total energy, with Z being the charge number and Φ being the equilibrium electrostatic potential. P φ is defined as P φ = MR 2φ − Zeψ. The last term associated with ν eff in equation (9) represents a simple Krook collision operator. The For energetic particles, we use a model that has an overall slowing down distribution in the particle energy space, combined with a Gaussian distribution in the particle pitch angle space [25,35] (10) Here, C(ψ) is the normalization function determined by the EPs density N h = dΓf 0 h , and corresponds to the crossover velocity of the EPs [35]. M i , M e , M h denote the mass of thermal ion, electron and energetic ion, respectively. ε h is the EPs birth energy. The electrons equilibrium temperature is denoted as T e . χ = υ /υ is the pitch angle, χ 0 is the center of the Gaussian function for the pitch angle dependence, C i are coefficients that can be chosen to specify various non-uniformity and anisotropy of the EPs distribution along the particle pitch angle [35].
One of the key parameters, which we shall carefully study in this work, is the width of the Gaussian function δχ. This parameter weakly depends on the particle energy (12) where δχ 0 is a parameter that we shall scan, and is often referred to as the width of the particle distribution in pitch angle, in this work. Theoretically, the full model (10) recovers the isotropic distribution model as δχ 0 → ∞. This has also been numerically demonstrated [35]. It is worthwhile to mention that MARS-K has two separate implementations, for the isotropic and anisotropic distributions.
The detailed description of the drift kinetic effects, implemented in MARS-K, can be found in [31]. The key physics element in this formulation is the mode-particle resonance operator. For thermal particles with ω d being the bounce-orbit-averaged precession drift frequency of particles, ω * N and ω * T the diamagnetic drift frequencies due to the density and temperature gradients, respectively. ω E is the E × B drift frequency. ε k = ε k /T is the particle energy normalized by the temperature T. m is the poloidal mode number, q the safety factor, and l the harmonic number in the bounce orbit expansion. ϑ = 0 for trapped particles, with ω b being the bounce frequency of trapped particles. ϑ = ±1 for passing particles, with ω b representing the transit frequency in this case. For energetic particles, the resonance operator is written as [36] Here, ε max = 3.52 MeV for fusion born α-particles. For NBI generated fast ion ε max is generally a function of the minor radius ε = ε/ε max . N h is the EPs density, and It is important to note that the mode eigenvalue ω enters into the resonance operator in a self-consistent manner, leading to a eigenvalue problem that is non-linear with respect to ω. In our work, we focus on TAE triggered by precession drift resonances of trapped EPs. The bounce resonance terms, though present in the MARS-K form ulation and in the code implementation, is not included in the present study.
The kinetic contribution from trapped EPs is written as [37] where φ (t) = φ(t) − φ t is the is the periodic part of the particle phase variation along the toroidal angle, and · denotes the average over particle bounce period. Here we have introduced a scaling factor Ξ, that we shall scan between 0 < Ξ < 1. The kinetic effects of EPs are fully included at Ξ = 1. Scanning of Ξ is useful since: (i) it helps us to understand the transition from the fluid results to full kinetic results; and (ii) computationally this helps us to trace the proper eigenvalues in the complex plane. We point out that MARS-K solves the MHD-kinetic hybrid equations in a non-perturbative manner. Therefore, besides the usual Landau damping of thermal species and EPs, the model also accommodates the self-consistent coupling to other MHD spectra such as the Alfvén continua and the kinetic Alfvén waves (KAW). For KAW, the model includes the electron Landau damping, as well as the electron col lisions through a simple (Krook) collision model. However, for the large part of the investigations presented in this work, we do not include the kinetic contribution from thermal electrons into the nonperturbative computations. Therefore the radiative damping, due to coupling to KAW, is mostly absent. We also neglect the finite orbit width effect [38].

Equilibrium
We use the HL-2A toroidal geometry, which has a plasma cross-section being nearly circular. In figure 1(a), the shaded region represents the plasma, with the black contour showing the plasma boundary. The outer contour in figure 1 indicates the location of a (thin) resistive wall surrounding the plasma, at the minor radius b wall . It is well known that the radial location of a resistive wall has an important influence on the stability of the resistive wall mode-one kind of global instability [39][40][41][42][43][44]. Between the plasma boundary and the resistive wall is a vacuum gap. The whole space outside the wall is also assumed to be vacuum. The main parameters of the chosen equilibrium are: the major radius R = 1.65 m, the inverse aspect ratio = 0.207, the total magnetic field B = 1.27 T, the central electron density n e (0) = 1.8 × 10 19 m −3 , the electrons and ions have the same temperature, with the on-axis value of T e (0) = T i (0) = 1 keV. The safety factor has an on-axis value of q 0 = 1.8, and the edge value of q a = 7.8. Figure 1(b) plots basic equilibrium radial profiles, with the black curve showing the equilibrium pressure, normalized by the factor B 2 0 /µ 0 (µ 0 is the permeability of free space), the red curve showing the EPs pressure (subject to the same normalization as for the total equilibrium pressure), and the blue curve showing the safety factor q. In this study, we assume that the total equilibrium pressure is the sum of the thermal particles pressure and the EPs pressure. We also assume that the total pressure is fixed while varying the EPs kinetic contribution δW K by the scaling factor Ξ. This leads to the following definitions: where p h = p h /p th , and p eq = p h + p th the normalized plasma equilibrium pressure. β = 2µ 0 p/B 2 is the ratio of the plasma pressure to the magnetic field pressure. In the following calcul ations, we fix the EPs pressure β * = 0.61 while varying the fraction of EPs kinetic effect via the scaling factor Ξ.

Effect of free boundary condition on TAE instability.
Free boundary condition is generally adopted in this study. In figure 2, we compare the computed eigenvalues of the instabilities with that obtained by assuming the fixed boundary condition, for an equilibrium at high β N regime, where β N is the normalized β, defined as β N = βaB/I with a being the plasma minor radius, B the toroidal magnetic field, and I the total toroidal current. In this example, is the normalized β limit with an ideally conducting wall. We consider an anisotropic particle pitch angle distribution for EPs, with χ 0 = 0 and δχ 0 = 0.2 in equations (10) and (12). The EPs pressure is also fixed at β * = 0.61 in this section. The other main parameters are: the wall position b/a = 1.2, the wall eddy current decay time (for the n = 1 perturbation) τ w /τ A = 2 × 10 4 , the EPs birth energy ε h = 50 keV and density n h = 0.02n i . In our computations, the drift kinetic effect due to the precessional drift resonance of trapped EPs is taken into account. The kinetic contributions from trapped thermal particles are found to be negligible, due to the fact that the precessional drift frequency of trapped thermal particles is much smaller than the TAE frequency, resulting in very weak mode-particle resonances. The kinetic effect from thermal passing particles has little influence on the TAE instability in this study. The adiabatic contributions (non-resonant contributions) from both thermal and energetic particles are included. Figure 2(a) indicates that the trapped EPs have a strong destabilizing effect. With the free boundary condition, three branches of instability are triggered due to precessional drift resonances with trapped EPs (solid curves). The threshold values of Ξ c for triggering different branches are labelled as Ξ 1 c , Ξ 2 c , and Ξ 3 c respectively. Here, superscripts '1', '2', '3' denote the first, second and third branches, respectively. The first branch with fixed boundary condition is also plotted (dashed curve). For the first branch, the growth rate increases with the EPs kinetic effect, approaching a constant value at large Ξ. The corresponding real frequency of the mode, shown in figure 2(b), is nearly a constant between Ξ 1 c < Ξ < 0.5, and slightly increases with Ξ at Ξ > 0.5. This slight upwards shift of the mode real frequency at higher Ξ, for the most unstable branch, is probably due to the modification of the mode structure. Note that the MARS-K computed real frequency with  fixed boundary is slightly higher than that with free boundary, consistent with the finding from [45]. The more counterintuitive result is the comparison of the mode growth rate (for the first branch) between the two boundary conditions. Naively, one would expect that the mode is more damped by the fixed boundary, and thus requiring higher critical Ξ c (i.e. Ξ 1 c ) to excite the instability. However, the MARS-K computations show a smaller Ξ 1 c with fixed boundary. Systematic investigations with other parameter variations (e.g. scanning the ideal wall minor radius) or with other equilibria (e.g. a large aspect ratio plasma) consistently show the same trend. A close examination of the results shows a significant difference in the mode eigenfunction (to be shown in a later figure), between the free boundary and fixed boundary cases. In other words, different boundary conditions lead to the excitation of different eigenmodes by EPs. It is the difference in the eigenmode that makes it plausible that a less unstable instability is excited with the free boundary condition. This point is also demonstrated by an analytic toy model devised in appendix.
The other two branches have very small (γτ A ∼ 10 −4 ) mode growth rates in the whole Ξ range. The corresponding real frequencies of the modes are nearly constants versus Ξ. We examined the mode eigenfunction for these two branches and found that the associated plasma radial displacements are mainly localized near the plasma edge region, where the EPs pressure is very low as shown in figure 1(b). During the Ξ scan, the EPs pressure is fixed, explaining why the computed growth rate, as well as the real frequency, of these two modes is not very sensitive to the Ξ in the range of Ξ > Ξ c . Figure 3 shows the corresponding perturbed kinetic energy for the instabilities shown in figure 2. All the energy terms in this paper are normalized by the plasma inertia associated with the radial displacement alone. At the marginal stability point for the TAE, δW f = −Re(δW K ), where δW f is the positive MHD potential energy of the stable mode [46]. Since the fixed boundary condition enhances positive δW f (as compared with the free boundary condition), at the marginal point, in order to balance the relatively larger δW f for the fixed boundary condition, larger |Re(δW K )| is required, as demonstrated by figure 3(a).
Certainly, the conventional energy balance argument, based on the perturbative approach, cannot easily explain the above result. On the other hand, with the non-perturbative approach, we propose the following heuristic explanation. Our starting point is that, when the mode is initially stable without EPs, adding EPs tends to drive the mode; on the other hand, if the mode has freedom to change its own eigenstructure (i.e. nonperturbative approach), the mode will try to adapt itself to resist the triggering by EPs. Now, since the free boundary condition gives the mode more freedom, it will try to find a state that is more difficult for the EPs to trigger it. The change of the mode eigenfuction, between the fixed boundary and the freeboundary cases, is evident from the MARS-K computations.
In addition, since the TAE mode structure is sensitive to Ξ as shown below, the potential energy of the mode also varies during the Ξ scan. Hence, it is not trivial to use the standard TAE dispersion relation to judge the dependence of the mode eigenvalue on the real and imaginary parts of δW K . Figure 3(a) suggests that the discrepancy of Re(δW K ) between fixed and free boundary conditions becomes smaller with increasing Ξ. This can be explained as follows: (1) at larger Ξ, the mode structure in the core region is similar for two cases (to be shown in a later figure); (2) the EPs density mainly deposits in the core region s < 0.6 in this study, which implies that the kinetic contribution mainly comes from the core region based on equations (7) and (8). Consequently, Re(δW K ) is similar for the two boundary conditions at larger Ξ. Figure 3(b) shows that Im(δW K ) still has discrepancy between the two boundary conditions at larger Ξ. This may be due to that fact that the mode eigenvalues, at larger Ξ, do not exactly overlap as they may appear in the log-scale plot shown in figure 2(a). This difference in the mode eigenvalue will affect δW K via the resonance operator. These numerical results for the first branch are consistent with the theoretical framework of GFLDR [47,48]. On the other hand, for the other two branches, the perturbed kinetic energy is insensitive to Ξ. This should be due to the fact that the real frequency of these modes is mainly localized in the TAE gap, as will be shown later ( figure 5(a)). In what follows, however, we shall only focus on the first branch, since this is the most unstable branch. The other way of identifying the (stable) TAE branch, without resorting to energetic particle destabilization, is to study the mode response to high frequency external magnetic perturbations (TAE antenna). This is carried out using the MARS-F code [49], with results reported in figure 4. Here, an external active coil is placed at the minor radius of 1.3a, with a sensor coil being located at 1.2a. The plasma response, in terms of the perturbed magnetic field amplitude, is recorded at the sensor position while scanning the frequency of the ac cur rent flowing in the active coil. The field perturbation peaks at certain frequencies, matching the specific eigen-frequencies of (stable) plasma modes [50]. Here, we examine two of the largest peaks. The first peak appears at Ω c1 = 0.148, and the second occurs at Ω c2 = 0.151 (all frequencies are normalized by the toroidal on-axis Alfvén frequency). Compared with figure 2(b), we find that the frequency of first branch (at Ξ = 0.6) and the frequency of the third branch match well with Ω c1 and Ω c2 , respectively.
Yet another important way of identifying the TAE is to analyze the spectra of Alfvén continua and to compare the location of continuum gaps with the computed mode structure. Figure 5(a) plots the n = 1 Alfvén continua that calculated by the ideal MHD eigenvalue code GTAW [51] with the slow-sound approximation [52] for the case of β N = 2.299. The poloidal harmonics of the radial displacement for the first branch from figure 2, for the cases of Ξ = 0.3 and Ξ = 1.0, are plotted in figures 5(b) and (c), respectively. Here, the mode amplitude is normalized to unity. The free boundary condition is assumed. We emphasize that these unstable eigenfunctions are self-consistently computed by MARS-K, meaning that modification of the eigenmode structure, by kinetic effects from EPs, has been taken into account. For the first branch from figure 2, the real frequency of the mode lies in the TAE gap shown in figure 5(a). The other two branches from figure 2 also have their real frequencies falling into the TAE gap.
As Ξ approaches the critical value for triggering the TAE (Ξ = Ξ 1 c ≈ 0.3), the m = 2 and m = 3 components peak near q = 2.5, with the m = 3 being the dominant harmonic, as shown in figure 5(b). At Ξ = 1.0, however, the m = 2 radial displacement becomes dominant, with its peak location near the q = 2 surface. The intersection of the m = 2 and the m = 3 harmonics occurs near the q = 2.5 surface in this case, as shown in figure 5(c). Comparing figures 5(b) and (c), we find that at higher Ξ, the mode structure becomes more global for the m = 2 and m = 3 components. The radial location of peaks of the mode amplitude shifts inwards to the magnetic axis. This demonstrates that the TAE mode structure can be modified by the EPs effect. In addition, for both Ξ = 0.3 and Ξ = 1.0 cases, there is a finite radial displacement at the plasma edge, due to the free boundary assumption. In experiments, the TAE mode structure therefore may depend on the beam injection power. The transport level of EPs has a strong relation with the radial extension of the instability [54]. It is plausible that the TAE with broader mode structure, at higher Ξ, can induce higher transport level of EPs.    figure 5(b), is substantially different from the classical TAE structure as shown in figure 5(c). This agrees with the observation made in [54], which the TAE can have very different mode structure, but similar mode frequency. In experiments, accurate measurements of the TAE mode structure can be made by using electron-cyclotronemission radiometer (ECE) [55] and reflectometer [56]. Figure 6 compares MARS-K computed mode structuresthe radial displacement of the plasma-in the 2D space, for different boundary conditions at Ξ = 1.0. Both the real (figures 6(a) and (c)) and imaginary (figures 6(b) and (d)) parts of the displacement are plotted, although we point out that the real and imaginary parts are not uniquely defined-they can be traded off by a toroidal phase factor without changing the eigenfunction. Several observations can be made. Firstly, the mode has a clear ballooning characteristic. Secondly, comparison between figures 6(a), (b) and figures 6(c), (d) show that the mode is more global with the free boundary condition. In addition, in the low field side, the mode extends from the core to the plasma surface. The displacement is also clearly enhanced at the plasma boundary in high field side. This in turn may enhance the radial transport of fast particles near the plasma edge, and/or enhance the fast ion losses during the TAE bursting [14,19,20,[57][58][59][60][61][62][63][64][65][66][67]. Figure 7 shows the normal displacement associated with the n = 1 TAE, computed for the case of Ξ = 0.3. Free boundary condition is used. Comparison between figures 6(c), (d) and figure 7 reveals that the mode structures are clearly different for different Ξ. At Ξ = 0.3, the mode is mainly localized in the middle region of the plasma column. Since the width of the particle-wave resonance layer is proportional to the mode growth rate, the small growth rate at Ξ = 0.3 implies the narrower layer of the particle-mode resonance. In addition, the EPs kinetic effect also contributes to the modification of the mode structure. Finally, at Ξ = 0.3, the continuum damping may also have more contributions in inducing the singular behavior in the mode structure. All these factors lead to the generation of more localized mode structure as shown in figure 7.

Effect of equilibrium pressure on TAE instability.
It is well known that the TAE instability is strongly related to the pressure gradient of fast particles [53,54]. Previous study also shows that the appearance of TAE modes depends on the gradient of the total plasma pressure [68]. In the previous subsection, the TAE study is carried out with a relatively high plasma beta β N > β ideal-wall N . However, the TAE instability is sensitive to the range of β N [24,69,70]. Hence, it is worth studying the dependence of the TAE instability on the β N variation. In this subsection, we investigate the influence of equilibrium plasma pressure on the TAE instability.  figure 8(a). Overall, the TAE is unstable in the whole β N range in this study. Figure 8(b) shows that the real frequency of the mode generally decreases with β N , a result that is consistent with that obtained in [29]. Our results seem to be in-consistent with the conventional picture that the TAE is suppressed when β N exceeds a critical value.
We point out that the TAE we study here is different from the conventional TAE, in that the TAE here couples to the unstable ideal kink mode, since we operate on plasmas with β N exceeding the no-wall, even the ideal-wall, Troyon limits. Therefore, the ideal kink is intrinsically MHD-unstable in  most cases that we consider. The corresponding TAE mode thus has the growth rate which is comparable to that of the ideal kink mode (order of 10 −2 in the toroidal Alfvén unit). Therefore, one reason that the TAE growth rate increases with β N , is the pressure driven unstable ideal kink mode.
The other reason is associated with the change of the continuum damping with increasing β N . Figure 8(c) shows that, for the case of β N = 2.074, the mode real frequency is localized within the TAE gap, with vanishing continuum damping. However, for the case of β N = 1.44 (shown in figure 8(e)) the mode frequency intersects with the shear Alfvén continuum, resulting in continuum damping. Figure 8 suggests that the dependence of the mode growth rate on β N is mainly due to the change of the continuum damping.
The influence of the equilibrium plasma pressure on the n = 1 TAE structure is demonstrated in figure 9(a), where compared are three cases with β N = 2.299 (dotted curves), β N = 1.515 (dashed curves), and β N = 1.44 (solid curves), respectively. The peak perturbation of all three cases are normalized to unity. For the cases of β N = 1.515 and β N = 1.44   figure 9(b). However, for the m = 2 harmonic, the distance of mode inward shift is much larger than the shift of the q = 2 rational surface, which is nearly fixed at the same minor radius with varying β N .

Effect of energetic particle distribution on TAE instabil-
ity. In the results reported so far, we have assumed an anisotropic EPs equilibrium distribution in particle pitch angle space, with fixed parameters δχ 0 = 0.2 and χ 0 = 0. On the other hand, previous studies showed that TAE instability is sensitive to the distribution of the EPs in particle phase space [25,26]. In this subsection, we shall investigate how the variation of the width of the Gaussian distribution in the particle pitch angle space can affect the excited TAE instability. Here, we fix the equilibrium plasma pressure at β N = 2.299, the EPs pressure at β * = 0.61, and the fraction of the kinetic contribution of EPs at Ξ = 0.8. According to equation (12), the anisotropic model, at the limit of δχ 0 → ∞, recovers the isotropic model, as previously benchmarked in a study on the resistive wall mode [35].
We perform a similar Gaussian width scan here for the n = 1 TAE, with results reported in figure 10. We fix the center of the particle pitch angle at χ 0 = 0, corresponding to a normal neutral beam injection. With the anisotropic distribution (blue curves), the mode growth rate initially slightly increases with δχ 0 , then sharply decreases with increasing δχ 0 , after reaching a maximum at δχ 0 = 0.35. Furthermore, the mode growth rate stays almost a constant as δχ 0 > 1. The real frequency of the mode sharply decreases with increasing δχ 0 when δχ 0 < 0.57, and then slightly increases with δχ 0 , approaching a constant with further increasing δχ 0 . Basically, the eigenvalue of the mode remains nearly unchanged as soon as δχ 0 > 1. At the limit of δχ 0 = 0 and χ 0 = 0, all EPs are trapped with the same pitch angle. As δχ 0 increases, fraction of trapped particles decreases, with the corresponding increase of the passing particle fraction. At δχ 0 → ∞, the EPs distribution in particle pitch angle space is uniform. For quantitative comparison, we also run MARS-K using the isotropic model. The computed eigenvalue (black lines) of the mode agrees well with that of the anisotropic model at large Gaussian width (δχ 0 > 1), as expected. Figure 10 demonstrates that the growth rate of the TAE instability is sensitive to the width of the particle distribution in pitch angle. Since the fraction of trapped particles at the local radial position decreases with increasing δχ 0 in our model (with fixed δχ 0 = 0), the instability drive for the TAE is reduced at larger δχ 0 . As a result, the TAE growth rate decreases with increasing δχ 0 in the region 0.35 < δχ 0 < 0.57. The mode structure trends to have a singularity layer in the same region of the Gaussian width. Change of the mode structure can also self-consistently affect the kinetic contribution from EPs. For the chosen case, in the region δχ 0 > 1, the anisotropy of the EPs distribution does not have significant influence on the TAE instability. In other words, at δχ 0 > 1, the anisotropic model is effectively equivalent to the isotropic model, for the TAE studied in this work. Figure 11 shows the mode structure in the 2D R-Z space for the case of δχ 0 = 10 (effectively corresponds to isotropic pitch angle model). Compared with the plasma displacement at δχ 0 = 0.2 (shown in figures 6(c) and (d)), the mode structure becomes more localized for the case of δχ 0 = 10, and the 'singular' layer appears near the q = 2.5 flux surface. This shows that the eigenmode structure also depends on the EPs distribution in the particle pitch angle space. There are two possible reasons. Firstly, the mode real frequency falls into the continuum spectrum of the Aflvén waves during δχ 0 scan. The resonance with continuum waves induces the singular behavior at the resonance position. This resonant coupling to continuum waves also reduces the TAE instability. Another explanation, which is less likely but cannot be excluded, is the conversion of the mode to another type of discrete Alfvén eigenmode. A systematic study of the mode conversion is beyond the scope of the present work. Figure 12 shows the 2D distribution of the absolute value of the complex perturbed kinetic pressure in the R-Z plane, defined as | p| = Re 2 ( p + p ⊥ ) + Im 2 ( p + p ⊥ ). Two cases, (a) isotropic and (b) anisotropic with δχ 0 = 0.2, are compared. Here, the kinetic pressure includes both parallel and perpendicular components. With the isotropic particle pitch angle distribution ( figure 12(a)), the perturbed kinetic pressure is more localized at the certain radial position, which is determined by both the mode structure and the EPs distribution in real space, as expected from equation (7). With the anisotropic case ( figure 12(b)), the perturbed kinetic pressure has a broader distribution, reflecting the more global eigenmode structure as shown in figures 6(c) and (d). Meanwhile, for both cases, the perturbed pressure is largely distributed in the core region, with very small amplitude near the plasma edge. This is due to the assumed EPs distribution function (shown in figure 1(b)), which yields very small ∇f 0 j near the edge.

Effect of plasma resistivity on TAE instability.
We have so far neglected the effect of the plasma resistivity on the TAE instability. It is, however, well known that the TAE instability is determined by the balance between (destabilizing) drive and (stabilizing) damping. The damping can include the continuum damping [21,22,71,72] as well as the resistive dissipation [23,73,74]. The relationship between the mode damping rate and the plasma resistivity depends on the resistivity regime [23]. Basically, at low plasma resistivity, the damping rate is proportional to η. At relatively high resistivity, the damping rate is proportional to η 1/3 .
In this work, we find that the plasma resistivity can completely stabilize the TAE instability, at least when the mode is near the marginal stability point. Figure 13 plots the mode eigenvalue in the 2D parameter space (δχ 0 , 1/S), for a plasma with β N = 2.299. Here, S is the magnetic Lundquist number. The white dashed curve in figure 13(a) indicates the marginal stability. In the small δχ 0 region, where the TAE has relatively large growth rate, the plasma resistivity is found to have a negligible influence on the mode eigenvalue. In the region of δχ 0 > 0.6, the plasma resistivity starts to slightly stabilize the mode. The TAE instability is completely stabilized in the region of 0.8 < δχ 0 < 1, with sufficiently high plasma resistivity. The critical value of the plasma resistivity for suppressing the TAE instability decreases (i.e. increasing the Lundquist number S) with increasing δχ 0 , as evident from figure 13(a). The main reason of the disparity between our results and Kar conclusion [23] is the modification of the radial mode structure due to the non-perturbative approach adopted in our modeling. More specifically, a singular layer appears due to EPs modification with an isotropic distribution in the particle phase space. The width of the singular layer increases with the plasma resistivity, indicating enhanced dissipation by plasma resistivity. Kar work did not include kinetic effects, and hence no kinetic modification on the mode eigenstructure can occur.
The real frequency of the mode decreases with increasing δχ 0 , before reaching a local minimum at δχ 0 = 0.57, followed by a slight increase with δχ 0 , for the whole plasma resistivity regime as shown in figure 13(b). Basically, the TAE has relatively large real frequency in the small δχ 0 region, and remains nearly unchanged while varying the plasma resistivity. However, in the region of δχ 0 > 0.6, the real frequency of the TAE slightly decreases with increasing 1/S. These results, shown in figure 13, mean that the plasma resistivity play a important role in the EPs-destabilized TAE, and contribute a considerable damping in a certain δχ 0 region, for the case of an anisotropic EPs equilibrium distribution in particle pitch angle space.
The influence of the plasma resistivity on the n = 1 TAE structure is demonstrated in figure 14, where two cases with S = 1 × 10 8 (solid curve) and S = 2.5 × 10 5 (dashed curve) are compared, respectively. Comparing figure 5(c) (δχ 0 = 0.2) and figure 14 (δχ 0 = 1), we find that at higher δχ 0 , the mode structure becomes more localized and a singular layer in the poloidal mode structure appears near the q = 2.5 flux surface. Obviously, the width of singular layer in the mode structure increases with the plasma resistivity, as shown in subplots of figure 14. The width of singular layer increases, as the mode eigenvalue varies during the plasma resistivity scan. The notable difference of the TAE mode structure between δχ 0 = 0.2 and δχ 0 = 1 is due to the fact that the real frequency of the TAE at δχ 0 = 1 drops in the continuum spectrum, and the TAE resonance with continuum waves induces the singular layer at the resonance location. Hence, in the large δχ 0 region, as continuum damping effect  increases with plasma resistivity, the TAE becomes more stable, as shown in figure 13.

Effect of wall location on TAE instability. Previous results
show that compared with the case of the fixed boundary at the plasma edge, a free boundary has a stabilizing effect on the TAE instability, with assumption that the radial location of resistive wall is fixed. However, since it is well known that the radial location of a resistive wall has a significant influence on the global instabilities, such as external kink mode [39][40][41][42][43][44], we in this subsection shall investigate how the radial location variation of a resistive wall affects the TAE instability. Figures 15(a) and (b) show the computed mode growth rate and real frequency, respectively, for the TAE while varying the normalized radial location of resistive wall b wall . The anisotropic particle pitch angle distribution for EPs is assumed with χ 0 = 0 and δχ 0 = 0.2. The TAE growth rate initially increases slightly with b wall . Then with increasing b wall , the growth rate stays almost a constant after reaching maximum at b wall = 1.3, as shown in figure 15(a). Figure 15(b) shows that the real frequency of the TAE monotonically decreases with increasing the b wall , and approaches a constant when the wall is far away from the plasma surface. This result is consistent with the previous finding in [29]. Figure 15 indicates that the wall position has a minor influence on the mode eigenvalue. This influence is due to the fact that the resistive wall  can act as an ideal wall for the TAE instability. Since the TAE real frequency is too high, and the perturbed magnetic field is difficult to penetrate the resistive wall due to the skin effect. This feature is demonstrated through monitoring the variation of the mode eigenvalue with the wall decay time τ w .

Effect of EPs birth energy on TAE instability.
In this paper, the TAE instability is excited due to precessional drift resonance with trapped energetic ions. As a result, the TAE instability should depends on the birth energy of EPs, due to the fact that the precessional drift frequency of the EPs is proportional to the birth energy of EPs. On the other hand, the previous study also shows that the TAE instability, trigged by passing particles, depends on the beam energy [25]. However, in the previous part, all of the computations are carried out with the assumption that EPs birth energy is fixed. Hence, in this subsection, we shall investigate the influence of EPs birth energy on the TAE instability. Figure 16(a) plots the eigenvalue of the mode in the complex plane, while scanning the normalized birth energy C = ε h /T i of EPs from 1 to 120. Here, the thermal ion temper ature T i = 1 keV is fixed during the C scan. The free boundary condition is adopted. At C = 1, the mode is an unstable ideal kink modified by the EPs kinetic effect, which recovers the standard ideal kink mode at the limit of vanishing EPs kinetic contribution (i.e. Ξ = 0). The mode growth rate increases with the EPs birth energy, reaching its maximum at C = 15. When C > 15, the mode growth rate decreases with increasing EPs birth energy. In addition, when C < 17, the mode frequency increases with C. However, when C > 17, the mode frequency enters into the TAE gap and has a weak dependence on C. Here, we identify the mode with the frequency within the continuum gap as the TAE. The mode with frequency smaller than the lower boundary of the TAE gap, ω r τ A = 0.14, is labelled as the modified kink mode. Here, C = 17 is defined as the cross point between the TAE and the modified kink, in the EPs birth energy space. Figure 16(a) shows that the kink mode can convert to the TAE mode with increasing EPs birth energy. Figure 16(b) shows the perturbed energy corresponding to figure 16(a), with increasing EPs birth energy. It shows that, at the cross point C = 17, Re(δW K ) + δW f is almost equal to zero. The TAE with frequency above the accumulation point ω r τ A = 0.14 exists in the region Re(δW K ) + δW f > 0, which is consistent with the analytical prediction from the general fishbone like dispersion relation [46][47][48].   Figure 17 shows the 2D structure of the n = 1 ideal kink with C = 1. In this case, the real frequency of the mode is in the order of 10 −3 ω A , far away from the TAE gap, as shown in figure 16(a). Comparison between figures 6(c), (d) and figure 17 show that the mode structures are notably different for different EPs birth energy. At C = 1, the mode is mainly localized in the core region of the plasma column. The dominant m harmonics is m = 2, and the mode structure is a kind of kink modes which modified by the EPs. This result indicates that the mode structure also depends on the EPs birth energy, and the mode conversion from ideal kink to TAE occurs during the EPs birth energy scan. In experiment, the mode transition can be demonstrated by measuring the real frequency of the mode during scanning the voltage of the inject beam and simultaneously measuring the mode structure through soft x-ray, ECE [55] or reflectometer [56]. The previous study shows that the TAE instability triggered by passing particles is stabilized, when the birth energy of passing particles is less than a critical value [25]. However, in this study, there is no energy threshold for triggering the instability when EPs pressure is fixed. This is due to the fact that the TAE originates from unstable ideal kink in this work. Since higher β N drives more unstable ideal kink, thus also drives TAE.

Summary and discussion
In this work, we have performed non-perturbative MHDkinetic hybrid computations of the TAE instability and structure, utilizing the MARS-K code. Several important plasma parameters are scanned. In this study, the TAE is driven by trapped energetic particles generated by neutral beam injection. Our computational results reveal that a realistic free boundary condition, as opposed to the conventional assumption of the fixed boundary at the plasma surface, leads to a more global mode structure, with finite plasma displacement in the edge region. We find that larger EPs contribution (larger Ξ) also results in a more global TAE structure, if the pressure of EPs fixed. Interestingly, at the threshold Ξ value for the TAE excitation, the mode structure, with the m = 2 and the m = 3 dominant poloidal harmonics overlapping with each other at q = 2.5, is significantly different from that of the conventional TAE. The bulk plasma pressure induces a slight inward shift of the m = 2 poloidal harmonic. EPs with an isotropic distribution in the particle phase space induce a singular layer in the mode structure, and the plasma resistivity enhances the width of this layer. Furthermore, a mode conversion, from an modified ideal kink by the EPs kinetic effect to the TAE, is observed while scanning the birth energy of EPs. The mode structure of the ideal kink is more localized compared with that of the TAE.
The free boundary condition increases the critical value of the EPs kinetic effect to drive the TAE. A resistive wall slightly stabilizes the mode, when the wall is sufficiently close to the plasma surface. The bulk plasma pressure basically enhances the TAE instability and reduces the real frequency of the mode. The growth rate of the TAE is sensitive to the anisotropy (in the particle pitch angle space) of the EPs distribution. Generally, both the growth rate and the real frequency of the TAE increase with enhancing the anisotropy. Finally, near the marginal instability point of the EPs driven TAE, the plasma resistivity can fully suppress the mode, yielding a stable domain in the 2D parameter space of the plasma resistivity and the width of the anisotropic pitch angle distribution.
The conclusions are presented in this work maybe can provide some heuristic insights for ITER operation: (i) ITER has a free boundary and a vacuum gap between the plasma boundary and the first wall. Our findings thus can be useful in giving qualitative predictions for the TAE instabilities in ITER. (ii) Several auxiliary heating schemes, such as NBI and ICRH, are planned for the ITER operation. These heating scheme will produce significant amount of trapped anisotropic EPs. Our finding, that both the growth rate and the real frequency of the TAE increase with enhancing particle anisotropy, can provide useful suggestions on the adjustment the EPs distribution (e.g. by adjusting the NBI injection angle).
In the future, the non-perturbatively computed TAE structure can be compared with experimental measurements, for realistic tokamak equilibria. The effect of the TAE structure modification, by free boundary or by other factors as investigated in this work, on the fast ion redistribution should also be an interesting research topic. Moreover, MARS-K will be used to study Alfvénic chirping, avalanching and large abrupt events in one or another way. boundary condition, compared to the case of a free boundary. Consider a 1D eigenvalue problem γy = γ 0 ∂ 2 y ∂x 2 + T(y), (A.1) with γ 0 > 0 and T(y) being an excitation term. Boundary conditions (BC) are specified at x = 0 and x = 1, respectively, with fixed boundary at x = 0 : y(0) = 0. Consider either fixed ( y(1) = 0) or free ( y (1) = 0) BC at x = 1. Consider a simple excitation model such as T(y) = αγy, with α > 1. The general solution of the above equation is y(x) = a cos(δx) + b sin(δx). The BC of y(0) = 0 yields y(x) = b sin δx. The fixed BC at x = 1 leads to a solution of δ = δ 1 ≡ π, whilst the free BC at x = 1 leads to δ = δ 2 ≡ π/2. The corresponding eigenvalue is Without the excitation term (α = 0), the linear system is stable (γ < 0). With sufficiently large excitation (α > 1), the system becomes unstable. More importantly, the free boundary solution is less unstable than the fixed boundary solution. This is essentially because different eigenmodes are excited depending on the BC, even with the same triggering term.
One can also devise a model that is phenomenologically closer to the TAE driven by EPs, e.g.
with iω = γ , which can again give more unstable instability with the fixed BC.