Simulation study of interaction between energetic particles and magnetohydrodynamic modes in the JT-60SA inductive scenario with a flat q≈1 profile

Interactions between energetic particles (EPs), an internal kink mode, and other magnetohydrodynamic (MHD) instabilities in the inductive scenario of JT-60SA (scenario # 2) are simulated with MEGA, a global EP–MHD hybrid code. For this scenario, it was predicted by TOPICS, an integrated transport code that the internal kink mode can be unstable and the sawtooth relaxation results in a flat safety factor (q) profile with q ≈ 1 for r/a⩽0.6 . In this equilibrium, it is found in the simulation results that the stability of the internal kink mode depends strongly on the bulk plasma pressure gradient ( ∇Pb ). In the n = 1 simulations where the toroidal mode number is restricted to n = 0 and 1, the pressure-driven internal kink mode is dominant. In the presence of co-passing EPs generated by the negative-ion-based neutral beam (NB), these EPs transfer energy to the internal kink mode; however, the EP driving rate ( γh ) is much lower than the driving rate from the bulk plasma pressure gradient ( γP ). The mode’s frequency is less than 1 kHz because the toroidal orbit frequency (ω φ ) and poloidal orbit frequency (ω θ ) of the co-passing EPs are approximately equal within the q = 1 surfaces. This mode affects the EP and bulk plasma pressure redistributions. The feasibility of stabilizing the internal kink mode using trapped EPs is also investigated. It is found that the trapped EPs with energy 85 keV generated by the positive-ion-based NBs cannot stabilize the internal kink mode. Stabilization is observed when the injection energy is greater than 500 keV. In the multi-n simulations, where n⩽8 modes are retained, the most unstable modes are high n interchange modes with poloidal number m = n whose linear growth rates exceed that of the pressure-driven internal kink mode observed in the n = 1 simulations. The overlapping of these modes creates a stochastic magnetic field, leading to stronger EP and bulk plasma pressure redistributions than those observed in the n = 1 simulations. During the nonlinear phase, the transition from the high m = n modes to the low m = n modes is observed where the dominant mode is the m/n=1/1 mode with an internal kink-like structure. These low m = n modes are generated by the nonlinear coupling of the high m = n modes. The EP kinetic effect has a minor contribution to the dynamics of these nonlinearly generated m = n modes.

In the JT-60SA inductive scenario with the lower single null (scenario #2), the magnetic field strength at the magnetic axis, the averaged toroidal beta (β b0 ) and the targeted plasma current are 2.25 T, 6.5%, and 5.5 MA, respectively.The 10 MW negative-ion-based neutral beams (N-NBs), the 24 MW positive-ion-based neutral beams (P-NBs), and the 7 MW electron cyclotron heating are used in this scenario, in which the non-inductive current drive fraction is 0.5.TOPICS [10,11], an integrated transport code, predicts that the internal kink mode can be unstable because the safety factor at the magnetic axis is less than unity (q 0 < 1).This mode can cause significant energetic particle (EP) and bulk plasma transport, which can affect the neutral beam (NB) current drive and bootstrap current.The negative-ion-based neutral beam injection (N-NBI) system installed in JT-60SA consists of two off-axis N-NBs which are injected at different elevations.The first and second beams are closer and farther from the magnetic axis, respectively.The combination of these off-axis N-NBs and the sawtooth relaxation (Kamdomtsev's reconnection model [12]) results in a flat q ≈ 1 profile.This flat q ≈ 1 profile can affect the magnetohydrodynamic (MHD) stability and the interaction between EP and MHD wave [13,14].
EPs generated by a NB injection (NBI) and an ion cyclotron resonance frequency (ICRF) heating can interact with the internal kink mode, causing both stabilizing and destabilizing effects.The stabilization effect arises from the conservation of magnetic flux within the area enclosed by the toroidal revolution of the trapped EP's precessional drift [15][16][17][18][19][20][21][22].This conservation is ensured when the precessional drift frequency (ω D,h ) is much higher than the linear growth rate (γ) of the internal kink mode (ω D,h ≫ γ).This condition is only satisfied for high-energy trapped EPs.If ω D,h ∼ γ, the fishbone instability [14,18,21,23] can be destabilized through a resonant interaction between the m/n = 1/1 mode and the EP [16,[24][25][26].The effects of EPs on MHD modes are not considered in TOPICS.Since EP can have both stabilizing and destabilizing effects on the internal kink mode, it is important to investigate the kinetic contribution of N-NB and P-NB-generated EPs in JT-60SA.In this study, the interaction between these NB-generated EPs and MHD modes in the inductive scenario of JT-60SA will be investigated with MEGA, an EP-MHD hybrid code.MEGA code has been successfully applied to several tokamaks [27][28][29][30][31].This study will focus on the following aspects: (1) the interaction between EPs and the internal kink mode in an MHD equilibrium with a flat q ≈ 1 profile, (2) the excitation of other EP-driven MHD modes in JT-60SA, (3) the redistributions of EP and bulk plasma pressure (P h and P b ) by the internal kink mode and other MHD modes, and (4) the feasibility of stabilizing the internal kink mode with trapped EPs in JT-60SA.This paper is organized as follows.In section 2, the simulation model in this study is reviewed.In section 3, the MHD equilibrium and the initial EP distribution functions (f h0 ) are presented.In section 4, the n = 1 simulation results with the N-NB-generated EPs are discussed.This includes the linear properties of the internal kink mode and its interaction with EPs.In section 5, the feasibility of the stabilization by trapped EPs in JT-60SA is examined.In section 6, the multi-n simulation results will be presented.Lastly, section 7 is devoted to discussion and conclusion.

Simulation model
MEGA [27][28][29][30][31][32][33] is a nonlinear hybrid simulation code, where nonlinear MHD equations are used to describe the bulk plasma, while drift kinetic equations [34] are applied to EP. MEGA solves nonlinear MHD equations with dissipation terms, ) where ρ, ⃗ v, P, ω, ⃗ E, ⃗ B, µ 0 , γ, η, ν, and ν n are bulk plasma density, velocity, bulk plasma pressure, vorticity, electric field, magnetic field, vacuum magnetic permeability, adiabatic constant, resistivity, artificial viscosity, and diffusion coefficients, respectively.The subscript 'eq' represents equilibrium variables.These coefficients are utilized to maintain numerical stability by dissipating the small-scale structures into thermal MHD energy through equation ( 3).These field quantities are solved in the cylindrical coordinate; therefore, the issue regarding the singularity at the magnetic axis on the m/n = 1/1 mode can be avoided.
The contribution of the EP to the bulk plasma can be found in the third term on the right hand side (RHS) of the momentum equation (equation ( 2)) in terms of the EP current density ( ⃗ J h ).The subtraction of ⃗ J h from the total plasma current density ( ⃗ J) is equivalent to the bulk plasma current density.⃗ J h is obtained through the moment integration of the EP guiding center distribution function.The contribution from the E × B drift of the EP is neglected due to the quasineutrality condition.This model is valid if the EP density is negligible when compared to the bulk plasma density.The MHD equations are solved explicitly by the 4th order finite difference method for spatial differentiation, and the 4th Runge-Kutta method for time integration.The finite-Larmor-radius effects are included through four-point gyro-averaging.In this study, four-point gyro-averaging is considered only in the multi-n simulation (section 6) because the scale length of the fluctuation is smaller for high-n modes.It will not be used in the n = 1 simulation (sections 4 and 5).
The δf method is employed in this study.The equilibrium marker distribution function is initially distributed uniformly in phase space (R, ϕ, z, v 3 , v ∥ /v), where R, ϕ, and z are cylindrical coordinates, and v and v ∥ are particle velocity and velocity component parallel to the confining magnetic field, respectively.Uniform marker distribution in v 3 is used in order to increase the high-energy marker population and to equalize the phase space volume that each marker covers.Under this condition, the number of physical particles that are presented by a single marker is proportional to the product of the equilibrium distribution function and the Jacobian of phase space J = 2π R i where i is an index of marker particle.The distribution function is normalized by a factor α which is calculated through the EP density (equation ( 7)).

MHD equilibrium
In this study, scenario #2 of JT-60SA is considered.The plasma parameters are modeled during the flat top phase by TOPICS, an integrated transport code.In this code, the turbulent and neoclassical fluxes are given by the current diffusive ballooning mode model and matrix inversion, respectively.The time evolution of the volume-integrated plasma stored energy (W total ) calculated with TOPICS is shown in figure 1(a).W total is composed of bulk electrons, bulk ions, and carbon impurities.In the TOPICS calculation, carbon impurity is assumed to have the same temperature as bulk ions.The density of carbon impurity at the magnetic axis is 2.57 × 10 18 m −3 , while the density of bulk ions is 6.06 × 10 19 m −3 .The sawtooth relaxation can be observed with the time interval of 1 s.The MHD equilibrium used in this study is reconstructed at t = 95 s, which is the time slice before the sawtooth crash.The poloidal flux of the equilibrium magnetic field, and the total plasma beta (β 0 ) and the safety factor (q) profiles at t = 95 s are shown in figures 1(b) and (c), respectively.The total plasma beta shown in panel (c) is composed of bulk electron (β e ), bulk ion (β i ), carbon (β c ), and N-NB-generated EP pressure (β h ).
The EP pressure is defined by Anisotropic EP pressure was not considered in the MHD equilibrium calculation.The total plasma beta at the magnetic axis (β b0 ) and the average plasma beta (< β b >) are 10% and 6%, respectively.The safety factor profile is flattened within the range of 0.1 < r/a < 0.6, where the safety factor is approximately equal to 1 (q ≈ 1).At the magnetic axis, the safety factor (q 0 ) is 0.92.This flatness is due to the combination of the Kadomtsev's full reconnection model [8,35] and the off-axis N-NBIs.

Energy and pitch angle distribution functions.
In this study, we assume the energy and pitch angle distribution functions of the N-NB and P-NB generated EP in JT-60SA based on the neutral injection system.The N-NBI system consists of two co-tangential beams with different trajectories, one close to the magnetic axis and the other further away.The P-NBI system consists of two co-injected, two counterinjected, and eight near perpendicular-injected beams.Our focus is on the N-NB-generated EP.We did not solve the EP distribution function from the Fokker-Planck equation; instead, we assume an anisotropic slowing-down velocity distribution, expressed by the analytic function where τ sd , τ cx , v c , v h,inj , ∆v, Λ, Λ 0 , and ∆Λ represent the slowing-down time, the charge-exchange time, the critical velocity, the injection velocity, the standard deviation of the complementary error function (0.1v h,inj ), the pitch angle, the center of the pitch angle distribution, and the width of the pitch angle distribution function, respectively.For the slowingdown distribution mentioned above, the ratio of the slowingdown time to the charge-exchange time (τ sd /τ cx ) is set to zero, which means that the charge-exchange time is much longer than the slowing-down time.The pitch angle parameter (Λ) is defined by µB 0 /E h where µ, B 0 , and E h are a magnetic moment, the magnetic field strength at the magnetic axis, and the kinetic energy of the EP, respectively.We assume a value of Λ 0 = 0.206, the same as the Fokker-Planck calculation results shown in [29].The energy and pitch angle distribution function of the N-NB-generated EP is shown in figure 2(a).This N-NB energy and pitch angle distribution will be used in sections 4 and 6.
For the positive P-NB generated EP, we do not have information about its pitch angle distribution function.However, we are interested in the role of the trapped EPs, so we only consider the perpendicular P-NBs in this study.The center of the pitch angle distribution (Λ 0 ) is set to 1.1.For the velocity distribution, we use the slowing-down velocity distribution function (τ sd /τ cx = 0) as shown in figure 2(b) for realistic JT-60SA plasmas.The figure shows that the peak of the pitch angle distribution is slightly below 1.1, which is a result of the non-uniform spatial distribution of EPs.Specifically, the parameter Λ 0 is defined based on the magnetic field strength at the magnetic axis, and EPs with Λ 0 > 1 are more likely to be found in the low-field side.Since the on-axis spatial distribution of EPs is used for positive P-NB generated EPs, integrating over space results in a decrease in the EP distribution at Λ 0 = 1.1.
To conduct an energy scan of the stabilization of the internal kink mode with trapped EPs, we also consider the bump-ontail velocity distribution.The results obtained using the bumpon-tail velocity distribution will represent the upper limit of achievable stabilization.This can be achieved by adjusting the ratio of the slowing-down time to the charge-exchange time (τ sd /τ cx ) to have a finite positive value.In this simulation, we assume a value of 4 for this ratio.This value is sufficient to reduce the low-velocity EP population.In the trapped EP energy scan, in addition to the maximum trapped EP energy (E h,max ) of 85 keV, we also explore higher energy levels of 500 keV and 1 MeV, as a part of the feasibility studies of the internal kink stabilization by trapped EPs on JT-60SA.The EP energy and pitch angle distribution for these additional cases are shown in figures 2(c)-(e).These EP distribution functions will only be used in section 5.

Spatial distribution function.
The initial EP spatial distribution is based on the equilibrium isotropic EP pressure profile calculated with TOPICS.The isotropic EP pressure is equal to the difference between the total plasma pressure shown in figure 1(c) and the summation of the bulk electron, bulk ion, and carbon impurity pressures.This isotropic EP pressure profile is labeled as 'Off-Axis N-NB β h0 = 1.62%' in figure 3(b).This is our main calculation scenario.With the assumed EP energy and pitch angle distributions shown in figure 2(a), we can recalculate the parallel and perpendicular components of the EP pressure, P h,∥ and P h,⊥ .In this approach, we assume that the velocity-averaged EP parallel (T h,∥ ) and perpendicular (T h,⊥ ) temperatures of EP are constant, which means that the EP pressure and density profiles have the same shape.The P h,∥ /P h,⊥ ratio is roughly 6.45 based on the pitch angle distribution.The maximum amplitudes of P h,∥ and P h,⊥ at r/a = 0.5 are 1.82 × 10 −2 and 2.82 × 10 −3 , respectively.In addition to the off-axis N-NB case with β h0 = 1.62% (red solid line), we also consider four additional cases in this analysis.These additional cases are the MHD case (purple solid line), the off-axis N-NB case with β h0 = 0.79% (red dashed line with asterisks), the off-axis N-NB case with β h0 = 2.37% (red dashed line with filled circular markers), and the on-axis EP case with β h0 = 5.68% (yellow solid line).The two additional cases with off-axis N-NB are considered for the scan of the EP beta, while the onaxis EP case is included to study the effect of injection location.Due to the smaller plasma volume near the plasma center, the on-axis EP case has a higher EP pressure at its peak compared to the off-axis cases.The MHD case is included to investigate the MHD limit.Since the total (equilibrium) plasma pressure must remain constant, the bulk plasma pressure (P b ) is equal to the difference between the total plasma pressure and the EP pressure (P h ).When the EP pressure is included in the simulation, the bulk plasma pressure gradient near the magnetic axis is enhanced for the off-axis cases but flattened and reversed for the on-axis case.Figures 3(c)-(g) illustrate the ratios between P b and P total and between P h and P total , respectively.The red and purple regions represent 'P h ' and 'P b = P total − P h ', respectively.These panels emphasize that the total plasma pressure remains constant for all the cases.

Simulation results: n = 1 N-NB case
In this section, we present the n = 1 simulation results with the N-NB-generated EP where only the n = 0, including zonal flow, and n = 1 modes are considered.Any perturbations with n > 1 are filtered out.We consider the MHD case, the offaxis N-NB case with β h0 = 0.79%, the off-axis N-NB case with β h0 = 1.62%, the off-axis N-NB case with β h0 = 2.37%, and the on-axis EP case with β h0 = 5.68%.Since we are dealing with the n = 1 mode, the four-point gyroaveraging is neglected.The cylindrical grid resolution (r, ϕ, z) used in this section is (12 816 256).To analyze the mode structure, we implement flux coordinates.The field quantity is mapped form the cylindrical grid to the grid of the flux coordinate.We then perform a Fourier analysis poloidally and toroidally in the flux coordinate

Linear growth rate
The time evolution of the energy fluctuation of the n = 1 mode for each case is shown in figure 4. For simplicity, we begin our discussion with the MHD case.The linear growth rate (γω −1 A ) of the MHD case obtained from the linear fitting is 2.11 × 10 −2 .From the radial MHD velocity (v rad ) profile (figure 5(a)), the dominant mode is m/n = 1/1 and its v rad profile peaks at the magnetic axis.Additionally, the v rad profile is limited inside the q = 1 surface (r/a ≈ 0.6).This suggests that this mode is an internal kink mode.The time evolution of the n = 1 energy fluctuation for the off-axis and on-axis EP cases are also shown in figure 4, while the v rad profile for the off-axis and on-axis EP cases are shown in figures 5(b)-(d), respectively.The significant modification in γω −1  A are observed.As shown in figure 4, γω −1  A is enhanced in the off-axis N-NB cases but reduced in the on-axis EP case.To understand these changes in γω −1 A , we estimate the EP driving rate (γ h ), the driving rate from the bulk plasma pressure gradient (γ P ), and the driving rate from the magnetic energy (γ B ).These driving rates are estimated by performing a volume integration of the total EP energy transfer ( ⃗ J h • ⃗ E), thermal plasma fluctuation energy, and magnetic fluctuation energy, respectively.To convert these values into driving rates, we calculate the time derivative for each term and divide it by the total energy of the n = 1 mode.This method can provide a good estimate only when there is a single dominant mode because all other perturbations are included.The results of this analysis are shown in figure 6.In this figure, γ h , γ P , and γ B are represented with red, purple, and green horizontal bars, respectively.The estimated γω −1 A of each case is shown on the RHS of the y-axis.Differ from γω −1 A shown in figure 4, γω −1  A in figure 6. is equal to the summation of γ h , γ P , and γ B .These results are in good agreement with the calculation results obtained from the linear fitting shown in figure 4. From the figure, it is clear that the internal kink mode is driven unstable by the bulk plasma pressure gradient (γ P > 0) while the magnetic field has a stabilizing effect (γ B < 0).The enhancements of γω −1  A observed in the off-axis N-NB cases with β h0 = 1.62% and β h0 = 2.37% are due to the enhancement of the bulk plasma pressure gradient (γ P ).The co-passing EPs also contribute to the energy transfer because γ h > 0 is always positive; however, γ h is much smaller than γ P .For the off-axis N-NB cases with β h0 = 0.79%, γ P is lower than the MHD case.This contradicts with the enhancement of γω −1  A observed in figure 4.This inconsistency may be attributed to the limitation of this method since all perturbations are included.The summation of γ h , γ P , and γ B shown in figure 6 (0.0228) is lower than the value obtained from the curve fitting shown in figure 4 (0.0267).The reduction of γω −1  A observed in the on-axis EP case is due to the reduction of the bulk plasma pressure gradient.The EP driving rate (γ h ) is higher in the onaxis case compared to the off-axis cases because it has a higher ∇P h in the core region.

Mode profile
For the n = 1 v rad profile, both EP fluid and kinetic effects cause modifications to the profile.In the off-axis, N-NB case with β h0 = 1.62% (figure 5(b)), the mode profile is slightly broader than in the MHD case.This change is due to changes in the bulk plasma pressure gradient and the kinetic and fluid responses of EP.In the on-axis case (figures 5(c) and (d)), a significant profile modification is observed.At ω A t = 571 (panel (c)), the m/n = 1/1 mode profile shows two peaks near the magnetic axis and at r/a ≈ 0.4, with different frequencies of 1.32 kHz and −0.5 kHz, respectively.This indicates that the mode profile consists of two different eigenmodes.The linear growth rates (γω −1 A ) measured at these peaks are 1.42 × 10 −2 and 1.29 × 10 −2 , respectively.The difference in the sign of the frequency indicates that the modes at these two radial locations are driven by different energy sources.Positive and negative frequencies correspond to ion diamagnetic and electron diamagnetic directions, respectively.This implies that the negative radial gradient of EP pressure drives the mode observed near the magnetic axis.We can conclude from this information that the peak near the magnetic axis is a fishbone mode driven by passing energetic particles, while the mode in the outer region is an internal kink mode.The locations of these peaks do not reflect the shape of the eigenmode because these two modes are overlapped.At ω A t = 856 (figure 5(d)), the minimization of the radial MHD velocity amplitude in the core region results from the destructive interference between these two modes resulting in the off-axis peak shown in the figure.This is similar to the linear results shown in figure 4.8 of [36], where the resistive kink mode has a slightly higher linear growth rate than the precessional drift fishbone mode.

Mode frequency
In terms of the mode frequency, the frequency of the mode observed in this JT-60SA equilibrium with a flat q ≈ 1 profile is low.The maximum mode frequencies of the off-axis N-NB β h0 = 1.62% and on-axis EP cases are −0.48 kHz and 1.32 kHz, respectively.These values are much lower than those observed in [24][25][26].The cause of this low-frequency mode is attributed to the flat q ≈ 1 profile, where the toroidal and poloidal orbit frequencies (f ϕ and f θ ) of passing EPs confined within the flat q ≈ 1 surface are approximately equal (f ϕ ≈ v ∥ /R 0 and f θ ≈ f ϕ /q).

Kinetic analysis
For the kinetic analysis in this section, we focus only on our main scenario which is the off-axis N-NB case with β h0 = 1.62%.Figures 7(a The flat q ≈ 1 profile region is located on the RHS of the r/a = 0.6 line.In panel (a), the total energy transfer is negative, indicating that the EPs transfer energy to the m/n = 1/1 mode.Panel (b) shows EP redistribution, with the interface between positive and negative δf h suggesting that the interaction between EP and the internal kink mode causes EP redistribution in the horizontal direction.This is expected because the frequency of the mode is almost zero.It is known that the combination of particle energy and toroidal angular momentum (E ′ = E − ωP ϕ /n) is a conserved quantity in the presence of a single dominant mode [37].With f = −0.48kHz, the E ′ = Const.line has a near-zero slope.However, in this case, the slope is not but close to zero because a finite EP energy transfer is observed.Panel (c) shows the f ϕ and f θ of all simulated EPs during the linear growth phase, with the contour representing |δf h |.It is apparent that the resonant EPs are located along the f θ = −f ϕ line (orange dashed line).These results support our previous explanation for the low-frequency mode observed in our simulation.Lastly, the poloidal resonance number (L) which is defined by equation (10), is shown in panel (d) where the poloidal resonance number is −1 which is expected since f ≈ 0 and |f θ | ≈ |f ϕ |.The resonant frequency ratio is independent of energy because f θ ≈ |f ϕ |.As a result, the poloidal resonance number is equal to a toroidal mode number.This implies that most of the passing EP within the q ≈ 1 surface can resonate with the mode, and it also explains the energy transfer layer shown in figure 7(a).

EP redistribution
During the nonlinear phase, the internal kink mode can induce strong EP (P h ) and bulk plasma (P b ) pressure redistributions.The redistributed P h and P b caused by the internal kink mode observed in the off-axis case with β h0 = 1.62% are shown in figures 8(a) and (b), respectively.Both the pressure profiles are nearly flattened within the q ≈ 1 surface (r/a ⩽ 0.6).The flattening is weaker for P h , while stronger for P b .This is also supported by the fact that this mode is driven by ∇P b .These redistributions can affect non-inductive current drive, plasma rotation, and q profile control in JT-60SA.

Simulation results: trapped EP n = 1 case
The influence of trapped EPs on the stability of the internal kink mode in JT-60SA is considered in this section.In JT-60SA, trapped EPs can only be generated by the perpendicular injection of P-NBs (E inj = 85 keV).For the ideal MHD model, trapped EPs can stabilize the internal kink mode through the conservation of magnetic flux, but only if the EP toroidal processional drift frequency significantly exceeds the linear growth rate of the mode.In this section, we consider the scans of the EP beta (β h0 ) at the magnetic axis and the maximum EP energy (E h,max ).Since E h,max can affect the EP orbit and the EP pressure profile, we assume the EP pressure profile fixed for all E h,max values.The on-axis EP pressure profile is assumed, while the bulk plasma pressure (P b ) is assumed to be the same as the MHD case shown in figure 1(b).It is fixed for every β h0 and E h,max value because the linear growth rate (γω −1 A ) of the internal kink mode observed in this equilibrium depends strongly on the bulk plasma pressure gradient (see section 4.1).The center of the pitch angle distribution parameter (Λ 0 ) is set to 1.1 to maximize the trapped EP population.For E h,max , we consider 85 keV, 500 keV, and 1 MeV.We consider different ranges of β h,max for each E h,max value.
For the E h,max = 85 keV case with a slowing-down energy distribution, the simulation results indicate that the stabilization of the internal kink mode by trapped EPs does not occur as shown in figure 9.This is also true for the E h,max = 85 keV case with a bump-on-tail energy distribution.These results suggest that the precessional drift orbit frequency of trapped EPs with E h ⩽ 85 keV is insufficient.In contrast, the internal kink mode can be stabilized for the E h,max = 500 keV case with a critical β h0 of around 0.01 and for the E inj = 1 MeV case with a critical β h0 of around 0.011.To further validate these results, we calculated the averaged precessional drift frequencies of the trapped EPs with E ≈ E h,max for each case and plotted them in figure 9(b).We also plotted the linear growth rates of the MHD case (γω −1 A = 2.11 × 10 −2 ) and the off-axis N-NB case with β h0 = 1.62% (γω −1 A = 3.09 × 10 −2 ) in the same figure as a black horizontal dotted line and a black horizontal dashed line, respectively.Our calculation of the toroidal precessional drift frequencies is consistent with the stabilization results, as the stabilization is observed only when the toroidal precessional drift frequency exceeds the linear growth rate of the MHD case.When compares to the linear growth rate of the off-axis N-NB case with β h0 = 1.62%, the trapped EPs with E h,max = 1 MeV are still sufficient to stabilize the internal kink mode.In the actual experiment, the degree of stabilization may be lower because the trapped EP energy distribution is not the bump-on-tail distribution.However, it can be intermittently realized via the selective confinement caused by the sawtooth crash [38].

Simulation results: multi-n case
In this simulation, we consider toroidal mode numbers within the n ⩽ 8 range.We choose to focus on these modes because, during the linear and nonlinear growth phases, the most unstable modes are the n = 7 mode and the n = 1 mode, respectively.The simulation with n ⩽ 16 has also been conducted but no significant differences in terms of the dominant mode and the saturation time have been observed; therefore, it is excluded from this manuscript.In this section, we will analyze the MHD case, as well as the off-axis N-NB case with β h0 = 1.62% and the off-axis N-NB case with β h0 = 2.37%.The number of grid points in the toroidal direction increases from 16 to 128 to accurately capture the n = 8 mode.Since high m = n modes have a smaller scale length, the finite Larmor radius effect is considered.In the first part (section 6.1), we will discuss the simulation results during the linear growth phase.P h and P b redistributions will be discussed in section 6.2.In the third part (section 6.3), we will analyze the nonlinear evolution of the low m = n modes and their interactions with the EPs.

Linear phase
During the linear growth phase, the dominant modes observed in both the MHD case, the off-axis N-NB case with β h0 = 1.62%, and the off-axis N-NB case with β h0 = 2.37% are the high-n modes with m = n.The time evolution of the fluctuation energy for each toroidal mode number for these cases are shown in figures 10(a)-(c), respectively.In this figure, the fluctuation energy of the n = 0 component has a negative value.This is due to the flattening of the bulk plasma pressure gradient.Since the negative value cannot be plotted in the logarithmic scale, the −E 0 component is plotted in place.The most unstable mode for all three cases during the linear growth phase is the n = 7 mode.The most dominant poloidal component for each toroidal number is m = n mode.The linear growth rates of all m = n modes for all cases at ω A t = 170 are shown in figure 10(d).The linear growth rate of the m/n = 7/7 mode for the MHD case, the off-axis N-NB case with β h0 = 1.62%, and the off-axis N-NB case with β h0 = 2.37% are 6.6 × 10 −2 , 7.2 × 10 −2 , and 7.3 × 10 −2 , respectively.These values are much higher than the internal kink mode observed in the n = 1 simulation (γω −1 A = 3.09 × 10 −2 ).The presence of the finite EP pressure slightly enhances the linear growth rate of the m = n mode, likely due to changes in the bulk plasma pressure gradient and resonant interactions.According to panel (d), the linear growth rate increases with increasing toroidal mode number (n), reaching a maximum at n = 7.This is similar to the behavior of ballooning and interchange modes, in which higher n modes are more unstable.The slowing down of the linear growth rate increment at high-n may be due to the viscous and resistive terms in the MHD equation, which are second-order spatial derivatives that can affect high-n modes.
To better understand the nature of these modes, we plotted the radial MHD velocity (v rad ) profile of the n = 7 mode of the off-axis N-NB case with β h0 = 1.62% in figures 11(a) and (b).The m/n = 7/7 mode is located within the q ≈ 1 surface, where shear stabilization is minimized.This mode exhibits a weak poloidal coupling and is dominated by a single poloidal harmonic m = 7.In contrast to pure MHD modes, it also has a finite sine component that leads to a shearing profile.This  modification is due to the EP kinetic effects, suggesting that this mode is an interchange mode modified by the EPs.The total energy transfer between EPs and MHD modes is shown in figures 11(c) and (d).In panel (c), we calculated the total EP energy transfer during 64 ⩽ ω A t ⩽ 193 in µ space for each mode.The total energy transfer is negative, indicating that the EPs transfer energy to the MHD modes.The maximum energy transfer occurs between the EPs and the n = 7 mode.Panel (d) shows the total EP energy transfer between the EPs with µ = 6.69 × 10 −2 and the n = 7 mode in (P ϕ , E h ) space.We can see destabilizing and stabilizing regions at r/a > 0.5 and r/a < 0.5, which are represented by blue and red colors respectively.This may be due to the off-axis N-NB profile where the sign of the radial gradient of EP pressures reverses at r/a ≈ 0.5.Similar to the n = 1 simulation results, these resonant layers also exhibit the same properties due to the low frequency of the mode and f θ ≈ |f ϕ |.

EP redistribution
In the multi-n simulation, the redistributions of EP (P h ) and bulk plasma pressures (P b ) are more severe than in the n = 1 simulation.The EP pressure and bulk plasma pressure redistributions in the multi-n simulation are shown in figures 12(a) and (b) for the off-axis N-NB case with β h0 = 1.62%.Compared to the n = 1 simulation, the bulk plasma pressure is completely flattened within the q ≈ 1 surface.The flattening is weaker for the EP pressure.This strong redistribution is due to the overlapping of the multi-n modes, which leads to a stochastic magnetic field.The Poincaré plots of the magnetic field during the linear (ω A t = 32) and nonlinear (ω A t = 321) phases of the off-axis N-NB case with β h0 = 1.62% are shown in figures 13(a) and (b), respectively.Each individual magnetic field line has a unique color, with blue representing the field line initially traced at the magnetic axis.The transition from blue to red represents the radial outward shift of the initial traced location.In figure 13(a), the flux surfaces are well-defined.However, at the saturation point, the magnetic field lines become stochastic (as shown in figure 13(b)).The red magnetic field lines diffuse into the core region and cover a large portion of the plasma volume which leads to a strong transport from the core region to the edge region.The stochastic field last for several Alfvén time which is sufficient for the radial transport of EP.

Nonlinear phase
During the nonlinear phase of all cases, the m/n = 1/1 mode becomes dominant.Figures 14(a    shift toward the core region, but this inverse cascading transition is much more pronounced for the off-axis N-NB cases (figures 14(g)-(j) and (l)-(o)).
In the MHD case, the excitation of the m/n = 1/1 mode is attributed to the nonlinear coupling between high m = n modes.This is supported by figures 10(a) and (e), which shows that the slopes of the n = 1, 2, and 3 modes become steeper than the n ⩾ 4 modes during 200 ⩽ ω A t ⩽ 300.The nonlinear growth rates of the m/n = 1/1, 2/2, and 3/3 modes are 0.131, 0.130, and 0.118, respectively, which are approximately two times higher than the linear growth rate of the n ⩾ 4 modes (γω −1 A = 0.066).The radial profile of the low m = n modes initially coincides with the high m = n interchange modes at ω A t = 285 (figures 15(a)-(h)).Additionally, the profile of the low m = n modes exhibits a tearing parity that agrees with the findings of [39], which indicate that a nonlinear mode generated by linear modes with an interchange parity will also have a tearing parity.Additionally, both the safety factor and the bulk plasma pressure profiles are flattened during the nonlinear phase; therefore, there is no additional free energy from the bulk plasma pressure and current density gradients to drive these low m = n modes during the nonlinear phase.
The transition from high m = n modes to low m = n modes is more pronounced in the off-axis N-NB cases.These low m = n modes observed in the off-axis N-NB case with β h0 = 1.68% are also driven by MHD nonlinearity, as evidenced by the agreement between the nonlinear growth rate (figure 10(e)) and the radial position of the nonlinear mode prior to saturation (figures 15(i)-(p)), with that of the high m = n interchange mode.Unlike the MHD case, the EP pressure profile is not completely flattened after the saturation of the high m = n modes, indicating the presence of available free energy from the spatial, energy, and pitch angle gradients of the redistributed EP distribution function.To investigate the EP contribution, the radial MHD velocity profiles of the m/n = 1/1 mode and the total EP energy transfer during the nonlinear phase of the off-axis N-NB case with β h0 = 1.62% are calculated and shown in figure 16.In panel (a), the m/n = 1/1 mode has a single dominant poloidal harmonic and it is located in the region with q ≈ 1.In panels (b), the net energy transfer between EP and each toroidal mode number is calculated in µ space within 579 ⩽ ω A t ⩽ 836.According to figures 10(b), 14(g) and (h), this time interval is the time when the m/n = 1/1 mode reaches its maximum amplitude.We do not observe the net negative (destabilizing) EP energy transfer to the mode during this time interval.Instead, the net energy transfer is positive (stabilizing) for the 1 ⩽ n ⩽ 8 modes (see figure 16(b)).The net stabilization is small when compared to the saturation amplitude.The total EP energy transfer between the EPs with µ = 6.69 × 10 −2 and the n = 1 mode in (P ϕ , E h ) space within the range of 579 ⩽ ω A t ⩽ 836 is also shown in panel (c).No apparent net negative energy transfer is observed.These results imply that the destabilization of the m/n = 1/1 mode observed in the off-axis N-NB cases is also caused by the MHD nonlinearity.The reason that the profile of the m/n = 1/1 mode is more pronounced in the presence of the EPs is due to the enhancement of the linear growth rates and the saturation amplitude of the high m = n modes.Another auxiliary factor is the weak stabilization effect of the EPs which are higher for the high m = n modes.For the n = 0 mode, the net energy transfer is negative (destabilizing); however, the enhancement of the n = 0 mode such as a zonal flow by EP is not observed.This negative energy transfer is possibly due to the source terms introduced in equations ( 1)-( 6) to restore MHD equilibrium.

Discussion and conclusion
In summary, we have investigated the MHD stability and the interaction between EPs and MHD modes in the JT-60SA inductive scenario with a flat q ≈ 1 profile.We have found that the most unstable mode is not the m/n = 1/1 pressuredriven internal kink mode, but rather the high m = n interchange modes.The overlapping of these m = n interchange modes significantly redistributes the EP and bulk plasma pressure profiles.Similar to the flattening caused by the sawtooth reconnection model, the bulk plasma pressure profile is completely flattened within the flat q ≈ 1 surface; however, the EP pressure profile is not.This suggests the need for an EPreduced transport model in an integrated simulation code.In addition to the investigation of the EP and bulk plasma pressure redistributions, we also have pointed out several physical aspects of the MHD stability and the interaction between EPs and MHD modes in such a unique equilibrium.They are as follows: (a) n = 1 simulation: -The stability of the internal kink mode mainly depends on the bulk plasma pressure gradient.The linear growth rate of this mode increases (decreases) in the off-axis (onaxis) N-NB cases mainly due to the enhancement (reduction) of the bulk plasma pressure gradient in the core region.The resonant interaction between the EPs and the m/n = 1/1 internal kink mode is also observed.The EP drive is weak in the off-axis case but strong in the on-axis case.-The frequency of the mode is low because the EP toroidal and poloidal orbit frequencies are approximately equal.-Most of the passing EP within the q ≈ 1 surface can resonate with the low-frequency mode.
-The pressure-driven internal kink mode cannot be stabilized using the trapped EPs generated by the P-NB system in JT-60SA.However, stabilization is observed when the population of trapped EPs with an energy of 500 keV or higher is sufficiently high.This is consistent with the theoretical prediction that stabilization occurs when the precessional drift frequency of trapped EPs exceeds the linear growth rate of the instability.However, generating trapped EPs with such high energies using NBs is challenging due to the shine-through problem.Therefore, an alternative method such as ICRF heating can be a feasible candidate for generating the required high-energy trapped EPs.(b) Multi-n simulation: -The EP pressure profile is not completely flattened by the overlapping of the high m = n interchange modes during the nonlinear phase.This may be due to the effects of magnetic drift and finite Larmor radius.-The transition from high m = n interchange modes to the low m = n modes caused by the MHD nonlinearity is observed.The dominant mode is the m/n = 1/1 mode with an internal kink-like structure in the core region.In the presence of EPs, the nonlinearly driven m/n = 1/1 mode has a higher amplitude.This is due to the higher linear growth rate and the higher saturation amplitude of the high m = n modes.The EPs also have a net weak stabilization effect.This nonlinear transition from high m = n modes to low m = n modes has been observed in other experiments, such as the sawtooth crash in DIII-D [40] and the pressure collapse in large helical device (LHD) [41], which both occurred in discharges with near unity flat safety factor and rotational transform profiles in the core region.There are two main differences between our study and those in DIII-D and LHD: our safety factor profile is at unity, rather than slightly close to unity, and we have considered the effects of EPs in our simulation.
These simulation results suggest that the high m = n interchange mode significantly modifies the bulk plasma and EP profiles.The enhancement of shear stabilization is difficult due to the lack of on-axis N-NBI in JT-60SA.Electron cyclotron current drive (ECCD) may be a possible option, but the 110 GHz gyrotron's second harmonic's cutoff frequency is close to the operation density (7.5 × 10 19 m −3 ).Interestingly, the instabilities can be used as another controlling knob.The overlapping of these m = n interchange modes can increase the EP population in the core region, resulting in more uniform heating and current drive.

Figure 1 .
Figure 1.JT-60SA scenario #2 parameters: (a) the time evolution of the plasma stored energy (W total ), (b) the poloidal flux of the equilibrium magnetic field, and (c) the total plasma beta and safety factor profiles at t = 95 s.

Figure 4 .
Figure 4. Time evolution of the n = 1 fluctuation energy in the n = 1 simulation for the MHD, off-axis N-NB, and on-axis EP cases.The purple solid line represents the MHD case, the red dashed line with asterisks represents the off-axis N-NB case with β h0 = 0.79%, the red solid line represents the off-axis N-NB case with β h0 = 1.62%, the red dashed line with open circles represents the off-axis N-NB case with β h0 = 2.37%, and the yellow solid line represents the on-axis EP case with β h0 = 5.68%.

Figure 5 .
Figure 5. Radial MHD velocity profile of the n = 1 modes calculated in the n = 1 simulations during the linear growth phase for (a) the MHD case, (b) the off-axis N-NB case with β h0 = 1.62%, and (c) and (d) the on-axis EP case with β h0 = 5.68%.Panels (c) and (d) show the mode profiles at ω A t = 571 and 856, respectively, for the on-axis EP case.The mode frequency for each case is shown in the figure.

Figure 6 .
Figure 6.Comparison of the driving rates from EP (γ h , red bar), thermal energy (γ P , purple bar), and magnetic energy (γ B , green bar) for the MHD, off-axis N-NB, and on-axis EP cases in the n = 1 simulation.γω −1 A obtained from the summation of γ h , γ P , and γ B are shown on the right-hand side of the y-axis.

Figure 7 .
Figure 7. Interaction between EP and the internal kink mode observed in the off-axis N-NB case with β h0 = 1.62% in the n = 1 simulation: (a) total EP energy transfer during the linear growth phase (ω A t = 429) for EPs with µ = 6.69 × 10 −2 , (b) δf h during the nonlinear phase (ω A t = 1286) for EPs with µ = 6.69 × 10 −2 , (c) EP toroidal and poloidal orbit frequencies during the linear growth phase, and (d) poloidal resonance number (L).The black dotted lines in panels (a) and (b) represent the normalized radius (r/a).The contours in panels (c) and (d) represent |δf h | during the linear growth phase.The orange dashed line in panel (c) represents the f θ = −f ϕ line.
) and (b) show the total EP energy transfer and the perturbed EP distribution function (δf h ), respectively: panel (a) at ω A t = 429 (linear phase) and panel (b) at ω A t = 1286 (nonlinear phase) for EPs with µ = 6.69 × 10 −2 .The gray dotted lines represent the normalized radius (r/a).

Figure 8 .
Figure 8.Time evolution of the (a) EP pressure and (b) bulk plasma pressure profiles in the n = 1 simulation for the off-axis N-NB with β h0 = 1.62%.

Figure 9 .
Figure 9. (a) Dependence of the linear growth rate (γω −1 A , blue line) and the frequency (red line) of the m/n = 1/1 internal kink mode and fishbone mode on the maximum EP beta (β h0 ) and the maximum EP energy (E h,max ).(b) Dependence of the averaged precessional drift frequency of the trapped EP on E h,max and their comparison with the linear growth rates of the m/n = 1/1 modes of the MHD case (black dotted line) and the off-axis N-NB case with β h0 = 1.62% (black dashed line).

Figure 10 .
Figure 10.Time evolution of the fluctuation energy for each toroidal mode number calculated in the multi-n simulation.Panels (a)-(c) show the results for the MHD case, the off-axis N-NB case with β h0 = 1.62%, and the off-axis N-NB case with β h0 = 2.37%, respectively.E 7,max shown in panels (a)-(c) represents the saturation energy of the n = 7 mode.The energy fluctuation of the n = 0 component has a negative value; hence, −E 0 is plotted instead.Panels (d) and (e) compare the growth rate of each mode at ω A t = 170 and 230, respectively.

Figure 11 .
Figure 11.Radial MHD velocity profile of the n = 7 mode and the total energy transfer between MHD modes and EPs in the off-axis N-NB case with β h0 = 1.62% during the linear growth phase.Panel (a) displays the radial profile of the n = 7 mode for different poloidal harmonics in flux coordinates, while panel (b) shows the n = 7 mode on a poloidal plane.Panel (c) shows the total EP energy transfer in µ space for each toroidal mode number within the time range of 64 ⩽ ω A t ⩽ 193.Panel (d) shows the total energy transfer between EPs with µ = 6.69 × 10 −2 and the n = 7 mode in the (P ϕ , E h ) space within the same time range.
)-(e), (f )-(j) and (k)-(o) show the time evolution of the radial MHD velocity (v rad ) profiles of the m = n modes for the MHD case, the off-axis N-NB case with β h0 = 1.62%, and the off-axis N-NB case with β h0 = 2.37%, respectively.Panels (a), (f ) and (k) show the mode profiles during the linear growth phase for these three cases, while panels (b)-(e), (g)-(j), and (l)-(o) show the mode profiles during the nonlinear phase at different time slices.The transition from the high m = n modes to the low m = n

Figure 12 .
Figure 12.Time evolution of the (a) isotropic EP pressure and (b) bulk plasma pressure profiles of the off-axis N-NB case with β h0 = 1.62% observed in the multi-n simulation.

Figure 13 .
Figure 13.Poincaré plots of the magnetic field during the linear and nonlinear phases of the off-axis N-NB case with β h0 = 1.62% in the multi-n simulation are shown in panels (a) and (b), respectively.

Figure 15 .
Figure 15.Radial MHD velocity profiles for the m = n modes of (a)-(h) the MHD case at ω A t = 285 and (i)-(p) the off-axis N-NB case with β h0 = 1.68% at ω A t = 257.These times are right before the end of the linear growth phase.

Figure 16 .
Figure 16.Radial MHD velocity profile of the m/n = 1/1 mode and the total energy transfer between EPs and MHD modes for an off-axis N-NB case with β h0 = 1.62% during the nonlinear phase.Panel (a) displays the radial profile of the first mode in flux coordinates.Panel (b) shows the total EP energy transfer within the range of 579 ⩽ ω A t ⩽ 836 in µ space for each toroidal mode number.Panel (c) displays the total EP energy transfer between the EPs with µ = 6.69 × 10 −2 and the n = 1 mode in (P ϕ , E h ) space within the range of 579 ⩽ ω A t ⩽ 836.