Effects of the resonance modification by electron cyclotron current drive on the linear and nonlinear dynamics of energetic particle driven magnetohydrodynamics modes in Heliotron J

This study investigates the effects of electron cyclotron current drive (ECCD) on the linear and nonlinear dynamics of energetic-particle (EP)-driven magnetohydrodynamics (MHD) instability (mode) in Heliotron J (HJ), using MEGA, a nonlinear hybrid EP-MHD simulation code. The effects of ECCD are included through the modification of the equilibrium magnetic field. Five MHD equilibria with plasma currents ( Ip ) of −2.00, −1.00, 0.00, 0.50, and 1.00 kA are considered, where negative and positive Ip represent counter and co-ECCD, respectively. To account for the EP finite charge-exchange loss observed in HJ, a bump-on-tail velocity distribution is assumed for EP. The study finds that the stabilization trends in terms of the linear growth rate ( γωA−1 ) are qualitatively similar to the FAR3D linear simulation and experimental results. In addition to the enhancements of the local magnetic shear and continuum damping by ECCD as proposed in previous studies, the study suggests that changes in the linear EP and shear Alfvén wave resonance are another important factor, especially for the n/m=1/2 mode. The study shows that in the currentless equilibrium, the high-velocity toroidicity-induced resonance between the n/m=1/2 mode and the co-passing EP ( L=−1 ) has the highest EP energy transfer rate, but is located slightly below the NBI injection velocity. The co-ECCD and cntr-ECCD can shift this resonance into higher and lower velocity regions, respectively. This can lead to resonance removal for co-ECCD, while cntr-ECCD can initially increase the number of resonant EPs for a small value of shift. If the shift of the resonance by cntr-ECCD is sufficiently large, the study observes stabilization due to the reduction in the EP population in the low-velocity region. The study also finds that changes in resonance affect nonlinear dynamics such as EP transport and frequency chirping. The convective transport of the co- and counter-passing EPs can be controlled, which can only be achieved for energetic-particle-mode due to the asymmetry in the resonance velocity between co- and cntr-passing EPs. The changes in the frequency chirping from the downward chirping branch to the upward chirping branch can also be explained by the changes in the resonance velocity caused by ECCD.

In this study, we focus on the effect of ECCD on the linear and nonlinear dynamics of EP-driven MHD modes.The application of ECCD can alter (1) local magnetic shear, (2) shear Alfvén continua (SAC), and (3) EP-shear Alfvén wave (SAW) resonance through changes in the rotational transform profile.In Heliotron J (HJ) [12][13][14][15], a low magnetic shear quasi-isodynamic optimized helical axis stellarator/heliotron with a toroidal field period of 4, the stabilization of the n/m = 1/2 energetic-particle-mode (EPM) and the n/m = 2/4 global Alfvén eigenmode (GAE) has been observed in both co-and counter-(cntr)-ECCD cases [1][2][3][4].The main proposed stabilization mechanisms are the enhancement of the local magnetic shear and continuum damping.The experimental results have also been qualitatively reproduced with linear simulation by FAR3D [4], a reduced MHD equation coupled with the Landau Closure models.The ECCD stabilization effects on the nonlinear dynamics (e.g.frequency chirping), the EP transport, and the changes in the EP-SAW resonance have not yet been widely investigated.In TJ-II, the change from continuous Alfvénic mode to the chirping mode by electron cyclotron resonance heating (ECRH) has been reported [16].In addition, Heliotron J's vacuum magnetic field has a near-zero shear with the rotational transform profile ranging from 0.54 to 0.57.This implies that the n/m = 4/7 and n/m = 4/8 magnetic islands can be generated by co-ECCD and cntr-ECCD, respectively.These islands can have a relatively large radial width [17,18].
MEGA [19], a hybrid EP-MHD code, has been used to study the linear and nonlinear interactions between EPs and MHD modes in tokamak [20,21], LHD [22], Heliotron J [23,24] (HJ), and CFQS [25].In MEGA, the MHD plasma is represented by nonlinear MHD equations in cylindrical coordinates, while EP is represented by the drift kinetic equation.For HJ, this code was qualitatively validated with the experimentally observed n/m = 1/2 EPM and n/m = 2/4 GAE in the low-beta, currentless equilibrium of the lowbumpiness configuration.This study also reports the importance of the free boundary condition in the modeling EP-driven MHD modes, where the interaction between EP and MHD mode is severely underestimated when a fixed boundary condition is applied.Since MEGA solves nonlinear MHD equations in cylindrical coordinates, it can treat an MHD equilibrium with magnetic islands.The objectives of these studies are (1) to investigate the modification of EP-SAW resonance by ECCD and (2) to investigate the nonlinear dynamics of EP-driven MHD modes in the ECCD-stabilized plasma of Heliotron J.
This paper is organized as follows.In section 2, the simulation model in this study is presented.In section 3, the simulation parameters including the HJ parameters, the MHD equilibria, and the initial EP distribution functions (f h0 ) are presented.In section 4, the simulation results for each scanned MHD equilibria are presented including the linear growth rates, mode profiles, mode frequencies, and EP-SAW resonances.In section 5, the nonlinear dynamics of the EP-driven MHD modes are discussed and compared with the experimental results are discussed.Lastly, section 6 is devoted to summary.

Simulation model
MEGA is a comprehensive MHD hybrid simulation code [19,20,22,24] that combines nonlinear MHD equations ( 1)- (6) to describe the behavior of bulk plasma and drift kinetic equations [26] to describe the guiding center of EPs.
The impact of EPs on the bulk plasma is included on the right-hand side of the MHD momentum equation ( 2), represented by the EP current density ( ⃗ J h ).The EP current density ⃗ J h is obtained through the moment integration of the EP distribution function.δf method is utilized in this study for EP.

MHD equilibria
In this study, we consider five MHD equilibria with a net toroidal plasma current I p = −2.00,−1.00, 0.00, 0.50, and 1.00 kA.The negative and positive values represent cntr-ECCD and co-ECCD, respectively.These equilibria are based on the low-beta plasma of Heliotron J in the low-bumpiness configuration, which is preferred in ECCD stabilization experiments due to the minimized bootstrap current.The averaged magnetic field strength (B 0 ), bulk plasma density (n 0 ), and bulk plasmas beta (β b0 ) at the magnetic axis are 1.18 T, 1.06 × 10 19 m −3 , and 0.36%, respectively.The MHD equilibria are prepared using HINT [17,18], a 3-dimensional MHD equilibrium calculation code that uses a time-dependent relaxation method and does not rely on the nested flux surface assumption.The EP pressure is not included in this MHD equilibrium calculation.In the MHD calculation, the current density profile for ECCD is approximated using (1 − s) 10 , where s represents the normalized toroidal flux.This profile is selected in order to replicate the localization of the ECCDdriven current.Since the plasma current profile is assumed, the scanned total plasma current range (−2.00 kA ⩽ I p ⩽ 1.00 kA) is not perfectly matched with the experiment.In the experiment, the plasma current is scanned from −0.8 kA ⩽ I p ⩽ 0.6 kA and −1.00 kA ⩽ I p ⩽ 1.00 kA for the n/m = 1/2 and n/m = 2/4 modes, respectively [3].The scanned total plasma current is approximately the same for the co-ECCD MHD equilibrium but it is 2 times higher than that of the experiment for the cntr-ECCD equilibrium.The reason that we extend the scanned range from −1.00 kA to −2.00 kA is to make the change in the rotational transform profile approximately the same for both co-and cntr-ECCD equilibria.The reason that cntr-ECCD equilibria require higher plasma current is possibly due to the expansion of the plasma volume observed in the cntr-ECCD MHD equilibria (see figures 1(a)-(c)).Since the plasma current profile is fixed concerning the normalized toroidal flux in our MHD equilibrium calculation, the change in the rotational transform profile will be weakened as the plasma volume expands in real space.The Poincaré plots of the equilibrium magnetic field and the rotational transform profile for each equilibrium are shown in figures 1(a)-(f ), respectively.The rotational transform profiles plotted in panel (f ) are plotted from the magnetic axis located at R = 1.35 m to the LCFS with Z = 0.00 m of the poloidal cross-section shown in panels (a)-(e).Magnetic islands with n/m = 4/7 are observed in the I p = 0.50 and 1.00 kA equilibria and their radial width is relatively large due to the low magnetic shear of Heliotron J. Finite-width magnetic islands with n/m = 4/8 are not observed in the I p = −1.00 and −2.00 kA cases because the ι/2π = 0.50 rational surface is located in the core region.In terms of the plasma volume, the volume decreases in the co-ECCD case but increases in the cntr-ECCD case.This may be due to the change in the poloidal magnetic field strength.
The formation of the n/m = 4/7 magnetic islands and the differences in the plasma shape between each MHD equilibrium can affect the spatial profiles of bulk plasma density, bulk plasma pressure, and EP pressure in real space.Since the bulk plasma pressure has been included in the MHD equilibrium calculation, the bulk plasma density profile is assigned to be consistent with the bulk plasma pressure.Within the O-point of the n/m = 4/7 magnetic islands, flat bulk plasma pressure and density profiles are assumed.This flattening is due to the rapid motion of thermal particles traveling along the magnetic field lines.The bulk plasma beta and density profiles plotted from the magnetic axis located at R = 1.35 m to the LCFS with Z = 0.00 m of the poloidal cross-section shown in panels (a)-(e) are shown in figures 1(g) and (h), respectively.

Initial EP distribution function
In Heliotron J, EPs are produced through both co-and counter-tangential NBIs of hydrogen.Each beam is characterized by the injection energy and the power of 28 keV and less than 700 kW, respectively.In prior studies, the energy and pitch angle distribution of EP has been assessed using charge-exchange neutral particle analyzers [27], revealing a bump-on-tail distribution.There are no diagnostic systems capable of directly measuring the spatial distribution function in HJ.However, it can be indirectly inferred from the density fluctuation of the EP-driven modes detected through beam emission spectroscopy where the n/m = 1/2 EPM and n/m = 2/4 GAE are observed near the edge region.Based on these data, it is reasonable to assume a strong EP pressure gradient near the edge region in our study [2,3].
The initial EP pressure profile is a critical factor in this calculation, as the EP-driving rate (γ h ) depends strongly on the EP pressure and its radial gradient.To simplify the calculation, we assume the EP pressure profile to be fixed in the flux coordinate for each MHD equilibrium; however, the difference in the shape of the magnetic flux surface and the magnetic islands can alter the EP pressure gradient in real space.To investigate this sensitivity, we consider 3 different initial EP beta (pressure) profiles with the steepest EP beta gradient (r ∇β h ,max ) at r/a = 0.55, 0.65, and 0.725, as shown in figure 2(a).These values are selected because the EP-driven MHD modes observed in the experiment are located near the edge region.The initial EP beta profile in real space for each case is presented in figures 2(b)-(d), respectively.The EP beta at the magnetic axis (β h0 ) is set to a fixed value.In our previous study [24], we used a value of β h0 = 1.10 × 10 −3 .However, we discovered that the EP-driven MHD mode was entirely suppressed in the case of I p = 1.00 kA.Since our interest is on the impacts of ECCD on the EP-driven MHD mode, we have increased β h0 to 1.50 × 10 −3 for this study.For equilibria with I p values of 0.50 and 1.00 kA, we assume that the EP beta (pressure) profile has a finite gradient within the O-point of the n/m = 4/7 magnetic islands due to the substantial magnetic drift of the EP.To achieve this, we employ spline interpolation, as illustrated in figures 2(b)-(d).
For the velocity distribution, we assume the bump-ontail distribution function as given in equation ( 7) [24,28] to account for the finite charge exchange loss in HJ [27].For simplicity, we consider only the full energy component of NBI.In equation (7), the variables f hs , v c , Λ, τ sd , τ cx , and v h,inj represent the initial distribution of EPs, the critical velocity, the pitch parameter, the EP slowing-down time, the chargeexchange time, and the NBI injection velocity, respectively.The NBI injection velocity in Heliotron J (HJ) is approximately 0.43v a0 , which is equivalent to 28 keV.The pitch parameter (Λ) is defined as µ B B 0 /E h , where µ B , B 0 , and E h represent the magnetic moment, the magnetic field strength at the magnetic axis, and the kinetic energy of EP, respectively.In this study, we set the ratios τ sd /τ cx , Λ 0 , and ∆Λ to 3.50, 0.20, and 0.35, respectively.These values are chosen to account for the bump-on-tail velocity distribution, tangential beam injection, and the helical magnetic axis of Heliotron J. Notably, these are the same values used in our previous study where we qualitatively reproduced the n/m = 1/2 EPM and n/m = 2/4 GAE [24].The assumed EP velocity distribution is visualized in figure 3.In the calculation, length, and velocity are normalized by ρ L = mv A0 /eB 0 and v A 0, respectively.

Linear simulation results
In this section, we investigate the linear properties of the observed MHD instabilities in each MHD equilibrium, including the spatial profile, frequency, and linear growth rate of the dominant modes.In section 4.1, the simulation results for the r ∇β h ,max = 0.550, 0.650, and 0.725 will be discussed and qualitatively validated with the experimental results shown in figures 10-12 of [3].The linear properties of the validated case will be discussed and compared with their SAC in section 4.2.
Lastly, the change in the linear EP-SAW resonance caused by ECCD will be discussed in section 4.3.
In the ECCD stabilized experiment in HJ [3], the n/m = 1/2 and n/m = 2/4 modes are observed in both the currentless, co-ECCD, and cntr-ECCD discharges.In the experiment, the maximum magnetic fluctuation amplitude of the n/m = 1/2 mode is strongly reduced with co-ECCD (I p > 0.00 kA); however, it is initially enhanced with the cntr-ECCD (−0.20 kA ⩽ I p < 0.00 kA) prior to the strong stabilization (I p < −0.20 kA).For the n/m = 2/4 mode, it is initially enhanced with the co-ECCD (0.00 kA < I p ⩽ 0.50 kA) prior to the strong stabilization (I p > 0.50 kA).Similarly, the enhancement of the n/m = 2/4 mode is initially observed (−0.25 kA ⩽ I p < 0.00 kA) with cntr-ECCD prior to the strong stabilization (I p < −0.25 kA).
In terms of dominant modes and stabilization trends, the case with r ∇β h ,max = 0.725 is the most similar to the experimental observations.However, five major differences between our simulations and the experiment are observed.(1) The co-ECCD case shows destabilization of the n/m = 2/3(4) toroidal Alfvén eigenmode (TAE).( 2) The cntr-ECCD case   (2, 3) The frequency of the n/m = 4/7(8) TAE observed in our cntr-ECCD case is within the same range as that of the n/m = 2/4 GAE reported in the experiment.From our simulation, the destabilization of the n/m = 4/7(8) TAE is observed only within a small parameter range.(4) The n/m = 2/4 EPM observed in the currentless MHD equilibrium for r ∇β h ,max = 0.650 and 0.725 is close to the extremum of the n/m = 2/4 shear Alfvén continuum (see figure 6(a)).A new EP pressure profile with a non-zero (weak) pressure gradient in the core region may be necessary to destabilize this mode.(5) The weak stabilization of the n/m = 1/2 EPM for the cntr-ECCD case observed in our simulation may be attributed to our assumption of a fixed EP spatial distribution.The actual EP beta profile and its maximum beta value may not be fixed in the experiment due to variations in the magnetic field structure and plasma volume (see figure 1).Based on these explanations, we will focus on the case with r ∇β h ,max = 0.725 which is the most similar to the experimental observations.The n/m = 1/2 mode has been observed in all MHD equilibria, and it is identified as an EPM for the I p = −2.00,−1.00, 0.00, and 1.00 kA MHD equilibria, as shown in figures 5(a)-(c) and (e).For the I p = 0.50 kA MHD equilibrium, it is identified as a GAE because its radial location and frequency are at the extremum of the n/m = 1/2 SAC, as shown in figure 5(d).It should be noted that the nested flux surface is assumed for the SAC calculation; therefore, the SAC at r/a ≈ 0.50 can be affected by the magnetic island.
Regarding the changes in the continuum damping, the slope of the n/m = 1/2 SAC is enhanced with cntr-ECCD, but it is reduced in the I p = 0.50 kA case.This change in continuum damping contradicts the stabilization trend of the n/m = 1/2 mode shown in figure 5(a), where γω −1 A of the n/m = 1/2 mode is enhanced and suppressed in the −1.00 and 0.50 kA cases, respectively.Since the change in the continuum by ECCD alone cannot fully explain the calculated results, the change in the EP-SAW resonance needs to be incorporated to provide a consistent explanation, which will be discussed in section 4.3.

n/m
The n/m = 2/4 mode is only observed in the MHD equilibria with I p = 0.00, 0.50, and 1.00 kA.In the I p = 0.00 kA MHD equilibrium, the n/m = 2/4 mode is identified as an EPM, as shown in figure 6(a).On the other hand, the I p = 0.50 and 1.00 kA (co-ECCD) MHD equilibria are identified as a TAE because the mode has a finite n/m = 2/3 harmonic, and both harmonics are located within the n/m = 2/3(4) TAE gap, as shown in figures 6(b)-(e).
Regarding the n/m = 2/4 mode, it is initially enhanced in the co-ECCD direction which is similar to the experiment.This may be attributed to the minimization of the continuum damping within the n/m = 2/3(4) TAE gap, as shown in figures 6(b) and (c).If the plasma current increases to 1.00 kA in the co-direction, stabilization is observed.In addition to the enhancement of the continuum damping, the reduction in the EP pressure gradient is also another possible factor since this TAE shifts towards the edge region where the EP pressure gradient is weakened (figures 6(d) and (e)).

Linear EP-SAW resonance modification with ECCD
This section discusses the modification of the EP-SAW resonance by ECCD.The alteration of the EP-SAW resonance can result from the change in the rotational transform (ι/2π) profile caused by ECCD, affecting both the SAC and the poloidal orbit frequency (ω θ ) of the EP.figure 8. shows the EP energy transfer rate (wdE/dt) in (v h , |v ∥ /v|) space during the linear growth phase for each MHD equilibrium.The negative contour represents the energy transfer from EP to the MHD mode.To calculate the EP resonance velocity, the assumption of the constant ι/2π value used in our previous study ( [24]) is invalid because our MHD equilibria have a finite magnetic shear.Instead, we trace the orbit of a test EP initially injected at the high-field side (HFS) and the low-field side (LFS) of the mode location.In real space, the HFS and LFS injected locations correspond to R = 1.20 m and 1.46 m of the HJ corner section shown in figures 1(a)-(e), respectively.These radial locations correspond to the mode location where the strong EP-SAW interaction is expected.The resonance velocity of the EP initially injected at the HFS and LFS are represented by the dashed and solid lines shown in figure 8, respectively.The red and blue colors represent the resonant velocity curves for the n/m = 1/2 and 2/4 modes, respectively.The poloidal resonance numbers (L, equation ( 8)) are shown only for the high-velocity toroidicity-induced resonances.The marked resonance is also the first sideband resonance.
As the n/m = 1/2 mode is not a coupling-type gap mode, there is no restriction on the symmetry between the toroidicityinduced resonances of the co-passing and counter-passing EPs.In the currentless MHD equilibrium (I p = 0.00 kA), the most important toroidicity-induced resonance is L = −1 resonance between the n/m = 1/2 EPM and the co-passing EP (figure 8(e)).This is followed by the L = 3 resonance between the n/m = 1/2 EPM and the counter-passing EP, as shown in figure 8(f ).The L = −1 and L = 3 resonance velocities are approximately 0.35 ⩽ v h /v A0 ⩽ 0.40 and 0.23, respectively, for the EPs with |v ∥ /v| = 1.00.The locations of these resonances are consistent with the negative EP energy transfer rate (blue region) calculated using MEGA.Furthermore, the accuracy of our analysis is confirmed by the L = 5 resonance between the n/m = 2/4 EPM and the cntr-passing EP, as shown in figure 8(f ).This L = 5 resonance is located in the high-velocity region (0.30 ⩽ v h /v A0 ⩽ 0.35); however, it has a low EP energy transfer rate than the L = 3 resonance between the n/m = 1/2 EPM and the cntr-passing EP at v h /v A0 = 0.23 located in the lower velocity region.This is possibly due to the lower amplitude of the n/m = 2/4 EPM compared to the n/m = 1/2 EPM during the linear growth phase.

co-ECCD (Ip
For the I p > 0.00 kA (co-ECCD) MHD equilibria, the L = −1 resonance of the co-passing EP shifts to a higher velocity region, as depicted in figures 8(e), (g) and (i).This shift is significant enough to exceed the NBI injection velocity, and it is observed only for the n/m = 1/2 mode.This shift can explain the significant stabilization of this mode seen in the I p = 0.50 and 1.00 kA MHD equilibria.In contrast, the shift is smaller for the cntr-passing EP, with a L = 3 resonance as shown in figures 8(f ), (h) and (j)).One of the primary reasons is that they are located in the lower-velocity region.The low-velocity EP has a smaller orbit width; therefore, it does not transit the core region where the strong ι/2π modification is observed.
For the n/m = 2/4 mode with v inj = 0.43v A0 , both the L = −2 and L = 5 resonances shift slightly towards higher velocity regions, but they still remain lower than v inj = 0.43v A0 .The changes in resonance velocity are smaller for this mode because it is a TAE in the I p = 0.50 and 1.00 kA MHD equilibria; therefore, it will follow the ι/2π = 2/3.5 surface.When the L = −1 resonance is slightly shifted into the lower velocity region, it can initially provide a destabilization effect because the number of the resonant EP is enhanced (In the I p = 0.00 kA case, the majority of the L = −1 exceeds v inj ).This observation is supported by the expansion of the negative energy transfer rate region in the I p = −1.00kA case (figure 8(c)), compared to the I p = 0.00 kA case (figure 8(e)).However, if the shift is too large, the number of resonant EPs can be reduced due to the bump-on-tail velocity distribution of the EPs.To verify our explanation for the removal of the L = −1 resonance, we conduct an additional simulation with increased NBI injection velocity (v inj ) for MHD equilibria with I p values of −2.00, −1.00, 0.00, 0.50, and 1.00 kA.In these simulations, we increase v inj from 0.43v A0 (28.00 keV) to 0.60v A0 (52.25 keV).Since we are focusing only on the L = −1 resonance, we apply a numerical filter to the perturbed EP pressure and retain only the perturbation within the N f = +1 toroidal mode family.

cntr-ECCD (Ip
We show the linear growth rate of the n/m = 1/2 mode for the v inj = 0.60v A0 case in figure 9. as the red solid line with asterisks.We find that the linear growth rate (γω −1 A ) increases for the I p = 0.50 kA MHD equilibrium, which corresponds to the shift of the L = −1 resonance towards the high-velocity region.As long as the shifted L = −1 resonance velocity remains below 0.60v A0 , the resonant EP population increases due to the bump-on-tail velocity distribution.However, when we increase plasma current to I p = 1.00 kA, we observe a reduction in γω −1 A .This is because the L = −1 resonance is close to 0.60v A0 , which results in a reduction in the number of resonant EPs.
In the cntr-ECCD direction, the linear growth rate decreases.This reduction is due to the shift of the L = −1 resonance towards the low-velocity region, and it is completely suppressed for the I p = −2.00kA case.
We also estimate the changes in the EP-driving rate (γ h , yellow solid line with asterisks) and MHD dissipation rate (γ d , green solid line with asterisks).They are shown in figure 9.

The change in γ h ω −1
A with respect to the plasma current agrees with the explanation we provided above.The change in γ d with respect to the plasma current is relatively smaller than γ h .

EP spatial redistribution
The main objective of the stabilization of EP-driven MHD instabilities is to suppress EP transport, making it important to investigate EP redistribution.Typically, we expect that EP redistribution should weaken as the linear growth rate decreases; however, EP transport is a nonlinear phenomenon affected by various factors, such as the mode number, frequency, amplitude, radial location, and resonance.In this section, we calculate the redistributed EP beta (β h ) during the nonlinear phase in each equilibrium and display the results in figure 10.This redistribution reduces EP drive and leads to the saturation of the EP-driven MHD instabilities.
In figures 10(a)-(e), the initial β h,0 is represented by the black solid line, and the red solid line represents the nonlinear phase.In all MHD equilibria, β h has a hollow beta profile because most of the EPs are high-velocity passing particles (bump-on-tail velocity distribution).Our previous study in [24].showed that these high-velocity EPs have a wide orbit width that enables them to interact effectively with the EPM in the outer region and transit the core region simultaneously.The convective transport of these high-velocity EPs leads to a significant reduction in the EP beta (pressure) in the core region.
For the I p = −2.00kA, −1.00 kA, and 0.00 kA equilibria, the trend of the redistributed β h follows the trend of the linear growth rate and the saturation amplitude shown in figure 4. and figure 11, respectively.However, the observed redistribution in the I p = 1.00 kA case (co-ECCD) is comparable to that of the I p = 0.00 kA case, despite having a much lower linear growth rate and saturation amplitude.This discrepancy can be explained by the change in resonance.According to our kinetic analysis discussed in section 4.3, the high-velocity toroidicityinduced resonance (L = −1) between the n/m = 1/2 mode and the co-passing EP is around 0.32 ⩽ v h /v A0 ⩽ 0.60 for EP with |v ∥ /v| = 1.00, while the high-velocity toroidicityinduced resonance (L = 3) of the cntr-passing EP is around 0.22 ⩽ v h /v A0 ⩽ 0.30 for EP with |v ∥ /v| = 1.00.This suggests that the n/m = 1/2 mode cannot cause the transport of the high-velocity counter-passing EPs.In contrast, the n/m = 2/4 and 2/3 modes have high-velocity toroidicityinduced resonances for both the co-passing (L = −2) and counter-passing EPs (L = 5), which implies that the n/m = 2/4 and 2/3 modes can transport both the high-velocity copassing and counter-passing EPs.
To confirm our explanation, we have separately plotted the redistributed EP beta (∆β h ) of the co-passing (∆β h,co , dashed line) and counter-passing (∆β h,cn , dotted line) EPs in figures 10(b), (d), (f ), (h) and (j).We observe that in the I p = −2.00,−1.00, and 0.00 kA equilibria, the convective transport of high-velocity co-passing EPs from the core region causes the majority of the redistributed EP beta, while the redistribution of counter-passing EP is weaker.In contrast, in the I p = 0.50 kA and 1.00 kA equilibria, the redistribution of counter-passing EP is stronger and comparable to that of the co-passing EP.These results provide support for our explanation of the changes in the EP-SAW resonance.

Frequency chirping
The frequency chirping can also be affected by changes in resonance due to ECCD.In the Heliotron J experiment, frequency chirping has been experimentally reported, especially for the n/m = 1/2 EPM [2,30].We compare the time evolution of the frequency spectra of the n/m = 1/2, 2/4, 2/3, 4/7, and 4/8 modes for each MHD equilibrium in figure 12. Panels (a)-(e) are for I p = −2.00,−1.00, 0.00, 0.50, and 1.00 MHD equilibria, respectively.The chirping branches of the n/m = 1/2, 2/4, 2/3, 4/7, and 4/8 modes are mainly dominated by a downward frequency chirping branch.At the current stage of our study, only the downward frequency chirping is validated with the experiment for the n/m = 1/2 EPM [2].The upward frequency chirping branch is only observed for the n/m = 1/2 EPM in the I p = −2.00kA MHD equilibrium.The change in the frequency chirping can be attributed to the change in resonance.To verify this, we need to evaluate the change in resonance for the L = −1 resonance, which is the resonance between the n/m = 1/2 EPM and the high-velocity co-passing EP through the toroidicity-induced resonance.This resonance is significant because it has the highest energy transfer rate during the linear growth phase for the I p = −2.00 and −1.00 kA MHD equilibria.
Figure 13 shows the perturbed EP distribution function (δf h ) during the early and later nonlinear phases for the I p = −2.00 and −1.00 kA cases.Panels (a)-(d) and (e)-(h) are for the I p = −2.00kA and −1.00 kA cases, respectively.Panels (a), (b) and (e), (f ) are plotted at the starting point of the nonlinear phase, while panels (c), (d) and (g), (h) are plotted at the later time of the nonlinear phase.The red solid lines represent the toroidicity-induced resonance between the EP and the n/m = 1/2 mode during the linear growth phase.The red dashed and dotted lines represent the toroidicity-induced resonance of the upward (ω U ) and downward (ω D ) chirping branches, respectively.For the I p = −2.00kA case, where the upward chirping branch is dominant, the high-velocity toroidicity-induced resonance of the upward branch is still lower than the NBI injection velocity.This implies that the upward chirping branch of the L = −1 resonance can interact with a higher number of EP populations when the L = −1 resonance shifts to the lower velocity region.In contrast, the upward chirping branch of the L = −1 resonance exceeds the NBI injection velocity for the I p = −1.00kA case.

Conclusion and discussion
In this study, we examined the influence of the rotational transform profile modification by ECCD on the linear and nonlinear dynamics of EP-driven MHD instabilities in Heliotron J by employing the MEGA code.The effect of the core-localized ECCD is included through the changes in the MHD equilibrium (e.g.magnetic shear and the rotational transform profile).These findings indicate that the EP-driven MHD modes are stabilized in both co-ECCD and cntr-ECCD directions, qualitatively consistent with the experiment and previous FAR3D linear simulation outcomes.Besides the previously proposed mechanisms of enhanced local magnetic shear and continuum damping, we propose that the change in the EP-SAW resonance is another crucial factor.The kinetic analysis conducted in this study suggests that the alteration in the EP-SAW  resonance significantly influenced the linear growth rate of the mode and its nonlinear dynamics.They are summarized as follows: (i) Co-ECCD can eliminate the high-velocity toroidicityinduced resonance (L = −1) between co-passing EPs and the n/m = 1/2 mode.In the absence of plasma current (I p = 0.00 kA), this resonance occurs at a slightly lower velocity than the NBI injection velocity.However, the modifications of the poloidal orbit frequency of the highvelocity EP and the mode frequency shift this L = −1 resonance into the higher velocity region, where it exceeds the NBI injection velocity in the co-ECCD (I p > 0.00 kA) cases.(ii) The enhancement of the n/m = 1/2 mode observed in the I p = −1.00kA cntr-ECCD case is attributed to the shift of the L = −1 resonance towards the lower velocity region.A slight shift of this resonance can increase the number of resonance particles because part of the L = −1 resonance exceeds the NBI injection velocity in the I p = 0.00 kA case.This mechanism can also explain the initial enhancement of the n/m = 1/2 mode observed in the experimental cntr-ECCD discharge.However, if the reduction in the L = −1 resonance velocity is too high, the number of resonance particles can be reduced because of the bump-on-tail velocity distribution.(iii) The effect of resonance shift on the n/m = 2/3(4) TAE observed in the I p = 0.50 and 1.00 kA co-ECCD cases is weaker, as the frequency of the mode is localized within the TAE gap.
The impact of ECCD on the nonlinear dynamics of EPdriven MHD instabilities is also significant.The changes in the linear EP-SAW resonance caused by ECCD can affect EP transport and frequency chirping dynamics.Normally, the reduction in the EP transport is expected when the linear growth rate (γω −1 A ) and the saturation amplitude (δv rad,max ) of the mode is reduced.Our results suggest that the change in the linear EP-SAW resonance caused by ECCD and the difference in the linear EP-SAW between different modes can also have a significant impact on EP pressure redistribution.This is supported by the comparable EP pressure redistribution observed in the I p = −2.00,0.00, and 1.00 kA cases.Lastly, the frequency chirping is also affected by the change in the linear EP-SAW resonance.The understanding of the change in EP-SAW resonances caused by ECCD may also apply to plasma diagnostics, especially in the case where the energy and pitch angle distribution is not experimentally available.Since SAC, EP-SAW resonance, mode number, and mode frequency can normally be obtained or calculated from the measured T e , T i , n profiles and Mirnov coils, the changes in linear growth rate and frequency chirping can qualitatively provide us the information about the EP energy and pitch angle distribution.However, the applicable operation range may be limited.
The assumptions we made when constructing the initial distribution function for the EP population may account for some of the discrepancies we observed between our simulation results and the experimental data.Specifically, our simulation did not account for the contributions of one-half and one-third of the NB energy components, nor did it consider how ECCD might alter the spatial distribution of EPs.These factors are known to influence both the linear and nonlinear (e.g.frequency chirping) behaviors of the mode we studied in the Heliotron J experiment.In terms of EP driving rate, the one-half and the one-third of energy components may increase the importance of the low-velocity EPs.The role of the first sideband resonance may be reduced.For the assumed Λ 0 , ∆Λ, and τ sd /τ cx values used in equation ( 7), these values can also be important to quantitatively reproduce the experiment.We have tried the case with a higher value of Λ 0 (Λ 0 = 0.4), the dominant mode becomes the n/m = 4/7 EPM (Have not been experimentally reported) and the n/m = 1/2 EPM is significantly weakened.The weakening of the n/m = 1/2 EPM can be easily understood from figure 8(e).As Λ 0 increases, EPs will populate at the lower value of |v ∥ /v|.At the lower value of |v ∥ /v|, the L = −1 resonance has a higher resonant velocity than the NB injection velocity.The scanned Λ 0 is excluded from this manuscript since it is out of the scope of this work.For ∆Λ, if a smaller value of ∆Λ is used, the linear growth rate will be more sensitive to the Λ 0 as the pitch angle distribution becomes more localized.For τ sd /τ cx , if a higher value is used, the role of the full-energy component will become more significant, and it will worsen our assumption where the one-half and one-third neutral beam energy components are neglected.If the small value of τ sd /τ cx is used, the velocity distribution will approach the slowing-down velocity distribution, and the n/m = 1/2 EPM is found to be weakened due to the reduction of the high-velocity EP.To quantitatively validate our simulation with the experiment, these factors must be carefully calculated.
The stabilization of the mode by the change in EP-SAW resonance observed in Heliotron J is only possible due to the bump-on-tail EP velocity distribution.However, in a burning plasma, where the maximum EP velocity exceeds the Alfvén velocity and the velocity distribution function is a slowingdown velocity distribution, the stabilization trend of the mode may differ, as demonstrated in the case where v inj = 0.60v A0 .These results highlight the need to extrapolate the proposed EP-driven MHD instability stabilization methods to burning plasma conditions.Additionally, the concept of asymmetry in the EP-SAW resonance between co-passing and cntr-passing EP for an EPM can be useful in plasma control with NBI.

Figure 1 .
Figure 1.Poincaré plots of the equilibrium magnetic field for the cases with the plasma current (Ip) of (a) −2.00,(b) −1.00, (c) 0.00, (d) 0.50, and (e) 1.00 kA.Panel (f ) displays the rotational transform profile of each equilibrium in real space, while panels (g) and (h) present the equilibrium bulk plasma beta and density profiles in real space.The profiles shown in panels (f )-(h) are plotted from the magnetic axis located at R = 1.35 m to the LCFS with Z = 0.00 m of the poloidal cross-section shown in panels (a)-(e).

Figure 2 .
Figure 2. Initial energetic particle (EP) beta profile (β h0 ).Panel (a) presents the initial beta profiles in the flux coordinate.β h0 profile of the r ∇β h ,max = 0.550, 0.650, and 0.725 cases are represented with filled circles, filled triangles, and filled squares, respectively.Panels (b)-(d) also show EP beta profiles for different equilibria in real space for the r ∇β h ,max = 0.550, 0.650, and 0.725 cases, respectively.The profiles shown in panels (b)-(d) are plotted from the magnetic axis (R ≈ 1.35 (m) to the edge (R = 1.55 (m) at Z = 0.00 m shown in figures 1(a)-(e).

Figure 5 .
Figure 5.Comparison of the radial magnetic fluctuation spectrum in the radius and frequency plane for the n/m = 1/2 mode observed in the r ∇β h ,max = 0.725 case during the linear growth phase with the N f = +1 shear Alfvén continua (SAC).Panels (a)-(e) are shown for the MHD equilibria with Ip = −2.00kA, −1.00 kA, 0.00 kA, 0.50 kA, and 1.00 kA, respectively.The red and grey solid lines represent the n/m = 1/2 shear Alfvén continuum and the other shear Alfvén continua within the N f = +1 toroidal mode family, respectively.The black contour represents the amplitude of the radial magnetic fluctuation.

Figure 6 .
Figure 6.Comparison of the radial magnetic fluctuation spectrum in the radius and frequency plane for the n/m = 2/4 and n/m = 2/3 mode observed in the r ∇β h ,max = 0.725 case during the linear growth phase with the N f = +2 shear Alfvén continua (SAC).Panels (a)-(e) are shown for the MHD equilibria with Ip = 0.00 kA, 0.50 kA, and 1.00 kA, respectively.The blue, violet, and gray solid lines represent the n/m = 2/4 shear Alfvén continuum, the n/m = 2/3 shear Alfvén continuum and the other shear Alfvén continua within the N f = +2 toroidal mode family, respectively.The black contour represents the amplitude of the radial magnetic fluctuation.

Figure 7 .
Figure 7.Comparison of the radial magnetic fluctuation spectrum in the radius and frequency plane for the (a) n/m = 4/7 and (b) n/m = 4/8 mode observed in the r ∇β h ,max = 0.725 case during the linear growth phase with the N f = +4 shear Alfvén continua (SAC).Both panels (a) and (b) are shown for the Ip = −1.00kA MHD equilibrium.The pink, purple, and gray solid lines represent the n/m = 4/7 shear Alfvén continuum, the n/m = 4/8 shear Alfvén continuum, and the other shear Alfvén continua within the N f = +4 toroidal mode family, respectively.The black contour represents the amplitude of the radial magnetic fluctuation.

4. 2 .
Mode spatial profile and frequency for the r ∇β h ,max = 0.725 case The comparison of the magnetic fluctuation spectrum of the n/m = 1/2 mode with the SAC within the N f = +1 toroidal mode family, the n/m = 2/4 and 2/3 modes with the SAC within the N f = +2 toroidal mode family, and the n/m = 4/7 and 4/8 modes with the SAC within the N f = +4 toroidal mode family for the r ∇β h ,max = 0.725 case are shown in figures 5-7, respectively.The solid red, blue, violet, pink, and purple lines shown in figures 5-7 represent the shear Alfvén continuum of the n/m = 1/2, 2/4, 2/3, 4/7, and 4/8 modes, respectively.These SAC are analyzed with STELLGAP[29].STELLGAP uses a nested flux surface equilibrium assumption; therefore, the effects of the magnetic island are not considered in the SAC of figures 5-7.

4. 2
.3.n/m = 4/7 and 4/8 modes.The n/m = 4/7(8) TAE mode is observed only in the I p = −1.00kA (cntr-ECCD) MHD equilibrium.The magnetic fluctuation amplitudes of these harmonics are comparable to each other within the TAE gap, as shown in figure7.Although this TAE has not been identified in experiments yet, its frequency is within the same range as that of the n/m = 2/4 mode.Further experiments are required to confirm the existence of this mode.

Figure 8 .
Figure 8.Comparison of the energy transfer rate (wdE/dt) between EP and SAW during the linear phase for each MHD equilibrium.Panels (a)-(j) are for Ip = -2.00,−1.00, 0.00, 0.50, and 1.00 kA, respectively.Panels (a), (c), (e), (g), (i) and (b), (d), (f ), (h), (j) are shown for co-passing and cntr-passing EPs, respectively.The solid and dashed lines represent the toroidicity-induced resonance velocity calculated by tracing the test EP orbit initially injected at R = 1.46 m (high field side) and R = 1.20 m (low field side), respectively.The red and blue colors represent the resonance of EP with the n/m = 1/2 mode and the n/m = 2/4 mode, respectively.The poloidal resonance numbers of the n/m = 1/2 and n/m = 2/4 modes (L) are shown only for the high-velocity toroidicity-induced resonances.

Figure 9 .
Figure 9.Comparison of the linear growth rate (γω −1 A , the red solid line with asterisks) of the n/m = 1/2 mode, EP driving rate (γ h ω −1 A , the yellow solid line with asterisks), and MHD dissipation rate (γ d ω −1A , the green solid line with asterisks) of the v inj = 0.60v A0 case.
< 0.00 kA) equilibria.In the cntr-ECCD MHD equilibria, the L = −1 and L = 3 resonances of the n/m = 1/2 mode shift slightly into the lower (figures 8(a), (c)) and higher (figures 8(b), (d)) velocity regions, respectively.The shifts of these resonances can have both stabilization and destabilization effects on the n/m = 1/2 EPM.To simplify the analysis, we focus on the L = −1 resonance (figures 8(a) and (c)) because its energy transfer rate is apparently larger than that of the L = 3 resonance (figures 8(b) and (d)).The effect of the shift of the L = −1 resonance on the stability of the n/m = 1/2 mode depends on the size of the shift.

Figure 12 .
Figure 12.Comparison of the time evolution of frequency spectra of the dominant modes caused by instabilities in each equilibrium.Panels (a)-(e) are for the Ip = −2.00,−1.00, 0.00, 0.50, and 1.00 kA MHD equilibria, respectively.

Figure 13 .
Figure 13.Perturbed EP distribution function (δf h ) in (v,|v ∥ /v|) space of the Ip = −2.00(a)-(d) and −1.00 (e)-(h) kA MHD equilibria.Panels (a), (c), (e), (g) and (b), (d), (f ), (h) show δf h of the co-passing and cntr-passing EPs, respectively.Panels (a), (b) and (e), (f ) show δf h at the starting time of the nonlinear phase, while panels (c), (d) and (g), (h) show δf h at the later time of the nonlinear phase.The red solid lines represent the toroidicity-induced resonance of the n/m = 1/2 mode during the linear growth phase.The red dashed and red dotted lines represent the toroidicity-induced resonance for the upward (ω U ) and downward (ω D ) chirping branches, respectively.