M3D-K simulations of beam-driven instabilities in an energetic particle dominant KSTAR discharge

We perform a systematic simulation study of energetic passing particle-driven instabilities in KSTAR using the kinetic-MHD hybrid code M3D-K. Linear simulation results show that the observed n = 1 mode in the early phase of the discharge is the low-frequency fishbone driven by energetic passing beam ions. The mode frequency computed is in a good agreement with the experimental measurement. Nonlinear simulations show that the frequency of the n = 1 mode jumps up to a higher value corresponding to the β-induced Alfvén eigenmode (BAE). In the later phase of the discharge, the simulated n = 5 mode is identified as a BAE in its linear phase. In the nonlinear phase, the n = 5 mode exhibits a similar frequency jump to a higher value of an energetic particle (EP) mode after mode saturation. Analysis of perturbed beam ion distributions in phase space shows that these new modes in nonlinear stages are driven by new resonances due to nonlinearly evolved beam ion distributions. Further simulations of a beam beta scan for the n = 5 mode show that the frequency jump disappears for a sufficiently small beam beta or beam ion drive. This result may explain the non-existence of frequency jump in the experiment. Finally, the impact of toroidal rotation on mode characteristics is investigated, showing that it has a marginal influence on EP driven modes.


Introduction
Energetic particles (EPs) from auxiliary heatings or fusion reaction in tokamaks would heat the background plasma through a collisional energy exchange.Especially the energetic alphas produced by the deuterium-tritium fusion Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.reaction are of crucial importance for a self-sustained burning in a fusion reactor.Meanwhile, they potentially destabilize shear Alfvén waves (SAWs) by inverse Landau damping via wave-particle resonant interaction [1][2][3][4].The characteristic speed of EPs in tokamaks is typically comparable to the Alfvén phase velocity in present tokamak experiments or future fusion reactors.This makes it possible for EPs to excite Alfvén eigenmodes (AE) or EP mode (EPM) by the radial gradient of an EP distribution via wave-particle resonances [5][6][7].In turn, these instabilities, with a characteristic spatial scale of EP gyroradius (mesoscale), would modify the EP distribution, enhance EP radial transport, and may even result in significant EP losses which is detrimental to the first wall of a fusion reactor [8,9].On the other hand, these instabilities could provide beneficial effects for helium ash removal from the core region or direct heating of the bulk plasmas (alpha channeling).Therefore, it is of great importance to investigate EP-driven instabilities and their effects on EP transport as well as on the stability and the confinement of background plasmas.
In a recent KSTAR experiment [10], low and high frequency chirping modes were observed with tangentially injected energetic beam ions.In order to better understand these beam-driven instabilities in the experiment, we have carried out a systematic simulation study of the beam ion-driven modes using parameters and profiles reconstructed from the KSTAR experiment.The global kinetic-MHD hybrid code M3D-K is used for this study [11].Simulation results show that, in the early phase of the KSTAR discharge, a low frequency fishbone is excited by passing beam ions with a frequency close to the measured one.Furthermore, in the later phase of the discharge, higher frequency β-induced AEs (BAEs) [12,13] are found to be driven unstable by beam ions in simulations.The calculated BAE frequencies are comparable to the experimental observations.Nonlinear simulations show that mode frequency jumps to a higher value for both the n = 1 fishbone and the n = 5 BAE.The details of the simulation results and a comparison with experimental observations as well as discussions will be given in later sections.
The fishbone instability was first observed in PDX [14].The mode has been named after the fishbone-like perturbed magnetic signal.The experimental discovery was followed by a theory development of the beam ion-driven fishbone.Fishbone driven by the spatial gradient of trapped EPs can be divided into two branches: the precessional drift branch (EPM branch) [15] and the diamagnetic drift branch [16].Chen et al [15] showed that the former branch comes from the internal kink mode destabilized by trapped EPs via waveparticle processional resonance.The frequency of the fishbone (EPM branch) is determined by the toroidal precessional drift frequency of trapped hot ions (defined as ω dh = −iv dh • ∇, v dh is the magnetic drift velocity of trapped hot ions) with the resonant condition ω ≈ ω dh .On the other hand, for the diamagnetic drift branch branch, the fishbone frequency is determined by the diamagnetic drift frequency ω * i of thermal ions.Later in 1986, the passing beam ion-driven fishbone was discovered in PBX [17].After its discovery, Betti [18] showed analytically that the finite orbit width effect is responsible for the passing beam ion-driven fishbone with the mode frequency determined by ω * i while the wave-particle resonance is given by ω = ω * i ≈ ω ϕ − ω θ .Subsequent works found that the EPM branch of the low frequency fishbone can also be excited by passing beam ions with the mode frequency determined by the same wave-particle resonant condition [19,20].Finally, the high frequency fishbone of the EPM branch was proposed to be excited by passing beam ions with a mode frequency determined by the wave-particle resonance of ω = ω ϕ [21].Most of the fusion-born alphas would be passing ones, so the research of energetic passing particle-driven instabilities is of crucial significance to the confinement of alphas in fusion devices.
Beta induced eigenmode was firstly observed in DIII-D with energetic ions from neutral beam injection (NBI) [12,13].The BAE frequency is comparable to that of the geodesic acoustic mode (GAM) and is inside the BAE gap of Alfvén continuum spectrum.The BAE gap is caused by the compressional response of the plasma to the propagation of SAWs in toroidal geometry [22,23].The finite compressibility is caused by the geodesic curvature of the equilibrium magnetic field in conjunction with the background plasma pressure.This produces a frequency upshift and leads to the BAE gap in the Alfvén continuum.Previous theory and simulations showed that the BAEs can be excited by both energetic ions [24][25][26][27][28][29] and energetic electrons [30].
This paper is organized as followed.Section 2 introduces the M3D-K code model and the major parameters and profiles used in simulations.Section 3 presents the simulation results in two phases of the KSTAR discharge.Section 4 presents a brief comparison of simulation results and experimental observations.Finally, in section 5, a summary is given.[11,31] which has been used to simulate EP driven instabilities [11,[32][33][34][35].It treats the plasma as a single fluid consisting of a set of full resistive MHD equations with additional dissipative terms.They comprise the continuity equation for the bulk plasma

Simulation model
the momentum balance equation the EP effects enter the equation (2) through the EP stress tensor P h , where P h = (P ∥ − P ⊥ )bb + P ⊥ I is the EP stress tensor in the Chew-Goldberger-Law form representing the EP effect [36], µ is the viscosity coefficient.And the adiabatic equation for the thermal pressure The thermal heat diffusion is not considered in this work.Equations ( 1)-( 3) are closed by Maxwell equations and the resistive Ohm's law Here, ρ denotes the bulk plasma mass density, v is the plasma velocity, P is the plasma pressure, J is the plasma current, E and B are electric and magnetic field, respectively.The δf method is used for reducingthe numerical noise with relatively small perturbations.The total distribution function is written as f = f 0 + δf = f 0 + wg, where f 0 and δf are the equilibrium and the perturbed distributions, respectively, g is the distribution function of simulation markers, and w = δf /g is the particle weight.The evolution of w can be derived from df dt = 0 and dg dt = 0 as The marker's orbits follow the guiding center equations in the drift kinetic limit where B 0 and δB are the magnitudes of the equilibrium and the perturbed magnetic fields, respectively, and b 0 is the unit vector in the direction of the equilibrium magnetic field.The angular bracket represents the gyro-average.B * and B * ∥ are given by and

distribution function of energetic particles
The slowing-down distribution function is assumed for the EP equilibrium distribution in our simulations [37].It is written as where erf(v) is the error function, v b is the beam injection speed, ∆v represents the width of the distribution near the injection speed, which is chosen to be ∆v/v b = 0.1, v c = (0.75 √ πm e /m i ) 1/3 (2T e /m e ) 1/2 is the critical slowing-down speed that is assumed to be a constant for simplicity, and c is a normalization factor.Ψ is the normalized poloidal flux defined as Ψ = Ψ −Ψ min Ψmax−Ψ min with Ψ min and Ψ max being the minimum and the maximum of Ψ, respectively, the angular brackets ⟨A⟩ represent an average of a physical quantity A over the particle orbit, and ∆Ψ represents the radial width of the distribution.Further, Λ = µB 0 /E is the pitch angle variable with Λ 0 = 0.2 the central pitch angle parameter and ∆Λ the width of the pitch angle distribution which is a function of particle energy E as given below [38]: where ∆Λ 0 = 0.25 is the width of the pitch angle distribution at the injection speed.Physically, equation (12) describes the broadening of a pitch angle distribution as beam ions slow down.Finally, since the toroidal angular momentum P Φ = eΨ + Mv ∥ RB ϕ /B is a conserved quantity, the orbit average of Ψ can be written as where e(M) is the beam ion charge (mass).We use the following approximate expression for passing particles.In our simulations, beam ions are mostly passing particles for tangential co-injection.

Parameters and profiles
We perform M3D-K simulations based on a KSTAR discharge #26759 displaying a variety of EP driven instabilities [10].
In this particular discharge, total 3.4 MW of NBI power was applied from t = 1.5 s, followed by 0.6 MW of electron cyclotron resonant heating at t = 2 s.We refer to [10] for details of discharge evolution.In this subsection, we describe main features of EP-driven modes in [10] for completeness.As reported in [10], the fast ion D a data indicate that strong passing beam ion density is formed in the core region.This significant fraction of EPs makes it suitable to investigate the onset and long-time-scale evolution of beam ion-driven instabilities in a rotating plasma.
The main parameters of the discharge [10] are as follows: the major radius R 0 = 1.8 m, the minor radius a = 0.45 m, the elongation κ = 1.7, and triangularity δ = 0.7, the toroidal magnetic field B 0 = 1.87 T, and the edge safety factor q 95 ≈ 4. The average injection energy of the deuterium beam is 90 keV.The corresponding normalized EP speed is v h /v A ≈ 0.45 with µ0ρ being the Alfvén speed.The normalized beam ion gyroradius is ρ h /a ≈ 0.07 at the injection speed.Figure 1(a) shows the spectrogram of the observed instabilities in the KSTAR discharge #26759.For convenience, the duration of the discharge is divided into three phases for analysis as indicated in figure 1: Phase A (2.2 s < t < 3.2 s), Phase B (3.5 s < t < 5.2 s) and Phase C (t > 5.3 s).In Phase A, multiple modes are excited with various frequencies.In particular, a lowfrequency chirping mode is clearly seen in Phase A. An expanded view of this mode is shown in the upper of figure 1(b).This  bursting chirping mode has m/n = 1/1 structure recognized from the T e fluctuation diagnosis by ECEI (electron cyclotron emission image).In addition, other modes with a frequency range between 100 ⩽ f ⩽ 200 kHz and mode numbers n = 3 and n = 4 are also observed.These modes gradually disappear at the end of Phase A. In Phase B, the low frequency chirping mode and other higher frequency modes mostly fade away.However, in Phase C, two high frequency chirping modes with toroidal mode numbers n = 1 and n = 5 appear as shown in the lower of figure 1(b).Fast frequency chirping of various AEs and EPMs have been widely observed in tokamak experiments, such as DIII-D [39], JET, JT-60U [40] and NSTX [41], and hole-clump pair generation in distribution function is found while the frequency is upshifting or downshifting [42,43].In this work, we focus on the analysis of the EP-driven low frequency chirping mode in Phase A and high frequency chirping modes in Phase C.
The experimental data are inputted to the TRANSP code to obtain the profiles shown in figure 2. It shows profiles of the safety factor (q), toroidal rotation velocity (v ϕ ), the total beta(β total ) and the beam ion beta (β h ) in Phases A and C. The x-axis (r/a) is defined as the square root of the normalized poloidal magnetic flux.We find the total central beta β total0 = 8.69%, central beam ion beta β h0 = 6.46%, and central toroidal rotation velocity v ϕ0 = 0.066 v A for Phase A, while β total0 = 5.54%, β h = 4.00%, and v ϕ0 = 0.045 v A for Phase C. Two notable differences between the two phases are the flattened q profile in the core region and a sharply decreased v ϕ in Phase C as compared to those of Phase A. Based on the given pressure and q-profiles and geometric parameters, we use the equilibrium code Variational Moments Equilibrium Code [44] to calculate a static equilibrium for each case as initial condition for M3D-K.Here we treat the two cases as two independent ones for analyses of beam iondriven instabilities.

Simulation results
We now present M3D-K simulation results of beam-driven modes based on parameters and profiles described in section 2. We will first describe results obtained without taking the plasma toroidal rotation into account in sections 3.1 for Phase A and 3.2 for Phase C. Effects of v ϕ are considered in section 3.3.In the M3D-K code, the extended MHD equations are solved as an initial value problem.The finite triangular element method is used for discretization in poloidal planes and finite difference method is used for discretization in toroidal direction.Typical numerical resolution parameters are 101 radial grid points and 400 poloidal grid points at the plasma edge.The number of poloidal grid points at each radial grid is proportional to radius.The time step length is typically chosen to be 0.01 τ A or smaller with τ A = R 0 /v A being Alfvén time.The number of simulation markers ranges from 12 millions to 40 millions depending on cases.The normalized viscosity parameter is set to be 1 × 10 −5 .A scan of viscosity shows that the growth rate decreases as the viscosity increases, and that the frequencies and mode structures vary little in both linear and nonlinear stages.We find these resolution parameters are sufficient for convergence of results.
In our simulations, real frequencies ω and growth rate γ are normalized to characteristic toroidal Alfvén frequency ω A , which is defined as A and 1000τ A ≈ 1.7 ms.

Simulations of phase A
Here we present simulation results for Phase A at the time slice t = 2.83 s.We use the nonlinear equations and do the numerical filtering at every time step by keeping all the n = 0 components fixed and keeping only a single toroidal mode number.We believe that this filtering process is equivalent to linear simulations although nonlinear equations are used for the MHD model.We use plasma profiles given in figure 2(a).We set the time step as 0.01 τ A , and use 12 millions simulation markers.Figures 3(a) and (b) show the evolution of the kinetic energy of the perturbed fluid motion (E k ) and the normalized frequency (ω/ω A ) of the n = 1 mode, respectively.The mode frequency is calculated using time evolution of a perturbation at a fixed point r = 0.3 with r being the normalized minor radius.In the linear stage, we observe that a low frequency mode appears with its mode structure shown in figures 4(a) and (c).The blue line in figure 3(c) indicates the mode frequency and the radial extent of the mode relative to the n = 1 Alfvén continuum spectrum (black curve).The Alfvén continuum is calculated using the NOVA code [45] (including the EP pressure effect).
It should be noted that the acoustic branches are excluded in this continuum calculation for clarity.Also the low frequency Alfvén branch is not included.However, the compressibility effect due to finite pressure is retained and as a result the shear Alfvén continuum is upshifted from zero at rational surfaces.This is possibly due to the resonance of the mode with the low frequency branch of Alfven continuum spectrum at a very small mode frequency and a weak growth rate.Future work will seek a more definitive mechanism for this particular mode structure.
In order to identify the main wave-particle resonance responsible for driving the low frequency n = 1 mode, we present figure 5 showing the color contour plot of the perturbed beam ion distribution in the phase space of (E, P ϕ ) with Λ = 0.2 in (a) the linear (t = 2800 τ A ) and (b) the nonlinear (t = 4200 τ A ) stages.The black line in figure 5(a) indicates the location of the wave particle resonance ω = ω ϕ + ω θ where ω is the mode frequency, ω θ and ω ϕ are the beam ion poloidal and toroidal transit frequencies, respectively.The general wave-particle resonance condition in toroidal geometry is give by where n is the toroidal mode number and p is an arbitrary integer.Since the perturbed distribution is localized near resonance lines for modes driven by the wave-particle resonant interaction, the result of figure 5(a) indicates that the low frequency mode is mainly driven by the p = 1 resonance.This result is consistent with an analytic theory of passing  particle-driven low frequency fishbone [19].In our M3D-K simulations, the negative sign of the mode frequency (see figure 5) corresponds to mode rotation in plasma current direction.The transit frequencies, ω θ and ω ϕ are positive and negative, respectively, for co-passing beam ions considered in this work.Near the q = 1 surface within which the mode is localized, ω θ and ω ϕ almost cancel each other in the wave particle resonance of ω = ω ϕ + ω θ .This is why this resonance drives a low frequency mode in accordance with the previous analytic theory [18,19].Finally, we note that there is a slight mismatch between the location of δf and resonant line at t = 2800 τ A (figure 5(a)).In general, the location of δf deviates a little from the wave-particle resonance line because δf also depends on other factors, such as the distribution gradient itself and the mode structure.Here the mode structure is mainly contained within the q = 1 surface.We believe that the mode structure is the main cause for the slight shift of δf location toward the center of the plasma.In summary, the linear simulation results show that the n = 1 mode is a low frequency fishbone driven by passing beam ions based on the analysis of mode structure and the main wave-particle resonance condition.
In the nonlinear stage, a careful observation of figure 3(a) shows that the mode first saturates at about t = 3400 τ A and decays for a short period before it grows again from t = 3700 τ A .During this period, the mode frequency jumps to a higher value of ω/ω A = 0.125.This upshifted frequency is located inside the BAE continuum gap around the q = 1 surface.The mode finally saturates for the second time at about t = 3900 τ A .The corresponding mode structure of this higher frequency mode is shown in figures 4(b) and (d).It is clear that the mode structure is similar to that of the low frequency mode found in the linear stage although their frequencies are distinct.From the mode frequency and mode structure, the higher frequency mode excited in the nonlinear stage (jumped mode) is likely a BAE mode.This is confirmed by the result of a parameter scan in the coefficient of specific heat Γ. Figure 6 shows that the frequency square of the jumped mode is linearly proportional to Γ which is consistent with the scaling of BAE frequency.Note that ω 2 BAE ≃ ω 2 GAM ∝ Γ.Instead of chirping observed in the experiment, jumping happening in the nonlinear stage is possibly due to the chosen EP distribution being not realistic and the lack of full damping physics in the code.Since the chirping physics is sensitive to the ratio of EP drive and damping as well as details of EP distribution function.
We now discuss the mechanism for the mode saturation and the frequency jump in the nonlinear stage.Figure 5(b) shows a color contour of the perturbed distribution function in the phase space around Λ = 0.2 in the nonlinear stage (t = 4200 τ A ).The solid lines represent p = 1 resonance (black curve), p = 0 resonance (blue curve), and p = 1/2 fractional resonance (red curve) of the mode, respectively.The corresponding radial profiles of the total distribution function are shown in figure 7 for particle energy around 60 keV at four different time slices: the start t = 0 (blue line), t = 3000 τ A before the first saturation (red line), t = 4000 τ A after the second saturation (yellow line) and t = 5000 τ A at the late saturation phase (purple line).The maximum amplitude location and the mode radial extent are indicated by the pentagram and a dashed black line, respectively.We observe that f is first flattened around P ϕ = −0.62 at t = 4000 τ A (yellow line) and the flattening region moves inward about P ϕ = −0.65 afterwards indicated by the purple line.The initial flattening of f causes the low frequency mode to saturate due to a reduced mode drive.The mode frequency jump and the excitation of the jumped mode, we believe, is driven by the p = 1/2 fractional resonance which induces the flattening of the distribution around P ϕ = −0.65.This is supported by figure 5(b) showing that the fractional resonance is aligned with the location of δf.It is known that the fractional resonance can be nonlinearly induced at large perturbations [46][47][48].Similar results of distribution flattening, frequency jump, and the fractional resonance are also observed in the simulation of the n = 3 mode in Phase A.
Finally, we study the sensitivity of the n = 1 result to the qprofile.Specifically, we consider the effect of q-profile shape and q 0 on the n = 1 mode within the experimental uncertainty of the profile.We consider two cases: the original q-profile (red line) and a new fitted one (blue line) as shown in figure 8(a).Note that the original profile is slightly reversed in the plasma core and the fitted one is monotonically increasing in radius.Figures 8(b) and (c) show the growth rate and the mode frequency, respectively, in the linear stage as a function of q 0 for these two profiles.We see that the growth rate is quite sensitive to q 0 .Furthermore, the trend of growth rate with increasing q 0 is opposite.This is due to the different shapes of the Figure 8. Two fitted q-profiles versus experiment data (circles).The growth rates (b) and frequencies (c) in the linear stages versus q 0 for two fitted q-profiles, respectively.q-profile.The mode frequency increases with q 0 when q 0 ⩽ 1, while it jumps to a higher value at q 0 = 1 for the fitted profile only.The mode structure of the low frequency mode in the linear stage becomes more centered in the core region as the q 0 value increase.However, simulation results show that the frequency jump in the nonlinear stage is insensitive to a small variation of a q-profile.The low frequency mode still jumps to a new one with a higher frequency of about 0.125 ω A for all the cases considered.This is understandable because the BAE frequency is insensitive to small variation of the q-profile in the plasma core.

Simulations of phase C
In this section, we consider beam ion-driven modes in Phase C of the KSTAR discharge.Recall that the n = 1 and n = 5 high frequency modes are observed with frequency down chirping.Therefore, we focus on a simulation study of these two modes.Time step sizes are 0.01 τ A and 0.002 τ A for n = 1 and the n = 5 simulations, respectively.Note that a smaller time step is necessary for the n = 5 simulation due to numerical stability.Figures 9(a) and (b) show the evolution of the perturbed kinetic energy of the n = 1 and 5 modes, respectively.The n = 1 mode has a much longer linear stage than that of the n = 5 mode due to a much lower linear growth rate (γ = 0.31% ω A for n = 1 and γ = 2.31% ω A for n = 5).The corresponding frequency evolution is shown in figure 9(c) for n = 1 and (d) for n = 5.We observe that, for n = 1, the mode saturates at about t = 4700 τ A and decays afterwards while the mode frequency stays nearly constant at about ω = 0.126.In contrast, the n = 5 mode saturates at an early time of about t = 700 τ A and decays a little afterwards while the mode frequency jumps from ω = 0.17 ω A to ω = 0.31 ω A at the mode saturation.Figures 10(a) and (b) show the mode frequencies and the radial extent of the modes relative to the continuum spectra for n = 1 and n = 5, respectively.Figure 11 shows the corresponding mode structures.We observe, for n = 1, the mode frequency is located in the BAE gap and the mode structure (figure 11(a)) indicates the dominant (1, 1) component.Thus the mode is identified as a n = 1 BAE.Likewise, the linear frequency of the n = 5 mode is also in the BAE gap around √ ψ = 0.5 while the linear mode structure (figure 11(c)) shows a dominant (m, n) = (6, 5) component.Based on the mode frequency and the mode structure study, the n = 5 mode is also identified as a n = 5 BAE.Furthermore, after nonlinear saturation, the n = 5 frequency jumps to a higher value which intersects with the Alfvén continuum at about √ ψ = 0.38 (red line in figure 10(b)).The corresponding mode structure of the jumped mode peaks around this intersection radius.Therefore, this higher frequency mode in the nonlinear stage is identified as an n = 5 EPM.Regarding the mode structure evolution shown in figure 11, it is interesting to note that the BAE-like mode in linear stage exhibits a pronounced shearing poloidal structure, while the EPM-like mode in nonlinear stage does not.This is because the shearing of the poloidal mode structure is mainly dependent on kinetic response of EP resonant drive [49].We believe that the shearing is closely related to mode growth rate due to EPs.In the linear stage, the growth rate is relatively large to form the shearing poloidal structure while in the nonlinear stage, the mode becomes saturated, and the shearing is much weakened due to a nearly zero growth rate.
We now discuss possible mechanisms for mode excitation and saturation as well as frequency evolution in the nonlinear phase of the two modes.Figures 12(a) and (b) show contour plots of the perturbed beam ion distribution and the main wave-particle resonance line in the phase space of (E, P ϕ ) at Λ = 0.4 for the n = 1 mode in the linear and the nonlinear stages, respectively.It is clear that the mode is driven by the p = 0 resonance in both linear and nonlinear stages.Note that there is no frequency jump in this case.Instead, there is only a slight upshift in frequency after saturation.It should also be noted that, for the n = 1 case, the main resonance occurs at a lower energy of about 20 keV and Λ = 0.4 based on the analysis of the perturbed distribution in the phase space.
However, for the n = 5 mode, situation is more complicated.Figures 12(c) and (d) show contour plots of the perturbed beam ion distribution and the main wave-particle resonance lines in the phase space at Λ = 0.2 for the n = 5 mode in the linear and the nonlinear stages, respectively.We observe that the p = 5 resonance is responsible for driving the mode  in the linear phase.In the nonlinear stage, the resonance is significantly broadened.Furthermore, a new p = 4 resonance (blue curve) appears.Lastly, the p = 4.5 fractional resonance also emerges.The dominant resonant energy range is 40-80 keV for the n = 5 mode.
We now consider the corresponding evolution of the beam ion distribution in the nonlinear stage of the n = 5 mode.Figure 13 shows that the distribution begins to be flattened around P ϕ = −0.57(red line) at t = 500 τ A , leading to the initial mode saturation around t = 500 τ A .Later, at t = 700 τ A when the mode saturates for the second time after a short period of the second growth, flattening of the distribution function becomes stronger with a much wider region.Meanwhile the second flattened region also appears at an inner radius of about P ϕ = −0.65.This new flat region is apparently caused by the new p = 4 resonance accompanied with a higher frequency mode.This suggests that the frequency jump is caused by new resonances driven by the gradient at inner radii.
The β h effects on the instabilities have been also investigated.We scan β h for the n = 5 while keeping the same thermal plasma beta (β thermal ). Figure 14 shows the normalized growth rate and the mode frequency as a function of β h .Naturally, the growth rate increases with β h , while the frequency and the mode structure nearly remain the same.The threshold β h for this linear instability can be obtained by extrapolation to zero growth rate, which is estimated to be 0.99%.We have also noticed the effect of β h on the frequency jump phenomenon in figure 15.With a decrease in β h , the frequency jump and the corresponding appearance of the jumped mode are delayed during the same simulation interval.The frequency jump even disappears when β h is considerably lower.This indicates that the low-frequency mode could be a BAE-like mode for such initial parameters and the jumped one is more affected by EPs, that is, an EPM-like mode.

Effects of toroidal rotation
In sections 3.1 and 3.2, the toroidal rotation has been excluded to focus on the identification of mode characteristics and driving mechanisms.In this subsection, we investigate the effect of finite v ϕ induced by tangential NBI.The M3D-K code is capable of retaining the toroidal rotation effect.
Figures 16(a) and (b) compare growth rates with and without v ϕ for several toroidal mode numbers in Phases A and C, respectively.The linear growth rate decreases for almost all cases except for n = 1 modes.For the n = 2 mode in Phase A with rotation, the mode is still stable as in the non-rotating case.This result is consistent with the experimental observation, showing the absence of the n = 2 mode.The simulated linear mode frequency increases with v ϕ for all the unstable modes while the mode structures are similar to those without rotation.
As for MHD modes, like kink and tearing instabilities, theories [50,51] have been proposed and simulation results have ascertained stabilization effects of rotation [52,53].For instance, deformation of the radial structure of toroidicityinduced AE from large v ϕ shear is deduced to be a reason for the stabilization effect by rotation in JT-60U [54] with counter-tangential NBI.However, no systematic and rigorous    Basically, it shows that it is MHD unstable with γ ≈ 0.6% at β h = 0 (Note that the finite frequency when β h = 0% arises from the toroidal rotation).This may explain the absence of  saturation of the kinetic energy since the MHD nonlinearity is neglected in the simulations.It was shown previous that MHD nonlinearity is important for saturation of a n = 1 mode driven by EPs [55].Therefore, no frequency jump appears for the n = 1 case in Phase A. Instead, the mode frequency settles down at ω/ω A ≈ 0.0628 ω A after a transient phase of frequency evolution.In Phase C (figures 17(c) and (d)), however, the mode does saturate and the frequency remains almost constant after nonlinear saturation with its value upshifted from the no-rotation value.
The results of the n = 5 case are shown in figure 19.At β h = 4% (figures 19(a) and (b)), we see that the mode frequency still jumps up to a higher value after saturation, which is similar to the case of no rotation.However, at a lower value  of β h = 2.3%, the frequency remains about the same in the nonlinear stage.This tendency is identical to the no-rotation case examined in figure 15.The higher frequency mode is not excited in the nonlinear stage during the same simulation period.
The f − P ϕ one-dimension distribution function approximately stays the same for the cases with relatively lower growth rates, such as n = 1, 2 in Phase C. In contrast, the n = 5 case in Phase C shows obvious flattening around the same radial location of the instability.A high linear growth rate results in the obvious reshaping of the distribution function around the mode radial location, which may induce new unstable modes driven by the enhanced gradient at both sides of the flattened location.
We have also carried out a scan of rotation speed to study the dependence of the linear growth rate on v ϕ for the n = 1 mode in Phase C. The rotation profile is fixed in this scan.Figure 20 shows that, as v ϕ increases, the growth rate firstly increases and then decreases.Mode structures and frequencies stay almost the same in this scan.We conclude that the rotational effects on the growth rate are modest within the range of experimental rotation speed for the n = 1 mode in Phase C.However, v ϕ has a significant stabilizing effect on n = 2, n = 3 and n = 4 modes in Phase C.This result is consistent with the experimental fact that beam-driven high frequency modes are not observed for n = 2, n = 3 and n = 4 in Phase C of the discharge.The physics origin for this intricate dependence of the linear mode characteristics on v ϕ is unclear.A theoretical understanding of this phenomenon will be a future subject.

Comparisons with experiment
We now compare our simulation results with experimental observations with respect to mode frequencies and their nonlinear evolution.Firstly, we consider a mode frequency comparison.The effects of v ϕ as well as diamagnetic frequency are included.Table 1 compares simulated frequencies of several modes with experimentally measured ones in Phases A and C. We include the effect of diamagnetic frequency ω * i = − nq eB0r ∂Ti ∂r with r being the radial location, T i the temperature of bulk ion at the q = 1 surface.The simulated frequencies of Phase A agree well with experimental measurements.For Phase C, the calculated frequencies of two modes are about 30% lower than the experimentally measured ones.Secondly, we compare frequency evolution.In Phase A, the frequency of the simulated n = 1 is nearly constant during nonlinear evolution.This is quite different from the observed frequency chirping in the experiment.The n = 1 mode frequency in Phase C also nearly remains constant while the experimental one exhibits a fast chirping.For the n = 5 mode in Phase C, the simulation results show frequency jump and chirping down after saturation at β h = 4%, but the frequency remains almost constant without jump and chirping at a lower beam beta value of β h = 2.3%.In summary, the simulated linear frequencies are similar to the experimental measurements.However, the simulated frequency evolution in the nonlinear stage is considerably different from the observation with respect to frequency chirping and jump.
The discrepancy in this frequency evolution may result from three sources.Firstly, the kinetic effects such as finite Larmor radius and thermal ion Landau damping are not included in our simulation model.Secondly, a simplified analytic beam ion distribution is used.Thirdly, the mode-mode coupling is not considered in our single-n simulations carried out in this work.These factors may have a considerable impact on mode stability and nonlinear behaviors of the simulated beam-driven modes and will be investigated in future.

Summary
Based on experimental data of an energetic ion dominant KSTAR discharge, we have carried out a systematic investigation of the energetic beam ion-driven instabilities using the kinetic-MHD hybrid code M3D-K.Linear simulation results show that, for an early phase of the discharge where EP pressure is relatively strong (Phase A), a strong n = 1 mode is excited resonantly by beam ions whose identity turns out to be the low frequency fishbone driven by energetic passing ions.The wave-particle resonance, mode frequency, and the mode structure computed from the code support this conclusion.The calculated mode frequency taking toroidal rotation into account agrees well with experimental measurement.Singlen nonlinear simulation results for the n = 1 mode show that the low frequency fishbone mode saturates due to a radial flattening of a beam ion distribution.Furthermore, after mode saturation, the mode frequency jumps to a higher value, converting its identity to a BAE.This transition is probably due to the linearly sub-dominant BAE becomes predominant with a newly evolved beam ion distribution in the nonlinear stage.At the late stage of the discharge where all plasma profiles are self-organized (Phase C), the linear simulation results show that both n = 1 and n = 5 modes observed in experiments are BAEs driven resonantly by beam ions.The calculated mode frequencies are similar to experimental values when v ϕ is accounted for.In the nonlinear stage, the simulation shows that the n = 1 mode saturates with the mode frequency nearly unchanged from the linear value.In contrast, the mode frequency of the n = 5 mode jumps to a higher value corresponding to that of an EPM after saturation.The EPM is driven by new resonances with a nonlinearly modified beam ion distribution.Furthermore, it is found that the frequency jump is sensitive to beam ion beta or beam ion drive.As β h decreases, the frequency jump is delayed until a sufficient deformation of the distribution function forms.For a sufficiently small β h , the frequency jump disappears and the mode frequency chirps down for the duration of simulation.
We have also investigated the effects of toroidal rotation systematically.It is found that the toroidal rotation generally stabilizes the beam ion-driven modes except for the n = 1 mode.Otherwise, toroidal rotation has little influences on mode structure or nonlinear evolution of the mode frequency.
Finally, regarding a comparison between simulation results and experimental observations, the mode frequencies for both cases generally show a good agreement.The main discrepancy lies in the nonlinear evolution of mode frequencies.Simulations show that the mode frequency jumps up to a higher value for the n = 1 mode without rotation in Phase A and for the n = 5 mode in Phase C.However, in the experiment, the mode frequency exhibits downward chirping for the n = 1 fishbone in Phase A and also for the n = 1 and n = 5 modes in Phase C.This discrepancy may come from differences between the analytic beam ion distribution and the realistic one in the experiment which can affect the beam ion drive and frequency chirping in the nonlinear stage.This possibility is supported by the results of a beam ion beta scan showing that the frequency jump is replaced by downward frequency chirping as β h decreases to a sufficiently small level.A more definite answer on this issue will be pursued in a future work.

Figure 1 .
Figure 1.(a) A spectrogram measured by Mirnov coils in the KSTAR discharge #26759.Expanded views of the n = 1 mode in Phase A and n = 1, 5 modes in the steady state (Phase (C) of the discharge are shown in (b).Adapated from [10].© IOP Publishing Ltd.All rights reserved.

Figure 2 .
Figure 2. Plasma profiles at two discharge phases.q is the safety factor, v ϕ is the toroidal rotation speed, β total (β h ) is the total (EP) pressure over the magnetic energy.Significant reduction in v ϕ and β h at Phase C can be seen in (b).

Figure 3 .
Figure 3.Time evolution of the perturbed kinetic energy (a) and the real frequency (b) of the n = 1 mode in Phase A. Time is normalization to the Alfvén time (τ A ). Fourier transformation amplitude is normalized to the maximum at each time.(c) The Alfvén continuum for the n = 1 mode (β h effect is included).Circles represent locations with maximum amplitudes and the colored lines show the radial ranges of the modes.Grey dashed line represents the safety factor of Phase A.
Figures 4(a) and (c) show the contour plots of the stream function U in poloidal cross-section and the radial structure of each poloidal component, respectively, in the linear stage.Here, U is the stream function of the incompressible part of the perturbed plasma velocity.The results show that the mode is consistent with the internal kink mode with a dominant (m, n) = (1, 1) component localized within the q = 1 surface.It should be noted that the (1, 1) mode shows a multi-ring like structure.

Figure 4 .
Figure 4. Poloidal ((a) and (b)) and radial ((c) and (d)) mode structures in linear ((a) and (c)) and nonlinear ((b) and (d)) stages.The radial location is defined as the square root of the normalized poloidal magnetic flux: r/a = √ ψp.|ϕmn| in (c) and (d) represents the amplitude of the mode with a certain m/n component.

Figure 5 .
Figure 5. Resonances in the linear (a) and the nonlinear (b) stages.Primary and secondary (half-integer) resonances are satisfied simultaneously in the nonlinear stage.

Figure 6 .
Figure 6.Adiabatic coefficient effect to the jumped mode.Γ is the adiabatic coefficient.

Figure 7 .
Figure 7. Evolution of the one dimension distribution function (f ) for the n = 1 mode in Phase A. The pentagram represents the location of P ϕ with the maximum mode amplitude and the dashed line denotes the mode range.

Figure 9 .
Figure 9. Evolution of the kinetic energy (a) and the normalized frequency (c) for the n = 1 mode.Those for the n = 5 mode are shown in (b) and (d).

Figure 10 .
Figure 10.Alfvén continua and radial mode extents showing that n = 1 (a) and n = 5 (b) modes are BAE modes in the linear stages (β h effect is included).Frequency jump in the nonlinear stage of the n = 5 mode changes the mode character from BAE to EPM.Circles show the radial locations with maximum amplitudes with lines representing the ranges.Grey dashed lines represent the safety factor of Phase A.

Figure 13 .
Figure 13.Evolution of the one dimension distribution function for the n = 5 mode in Phase C. The pentagram represents the location of P ϕ with the maximum amplitude and the dashed line denotes the mode range.

Figure 14 .
Figure 14.Normalized growth rate and frequency of the linear mode versus energetic particle beta (β h ).

figure 17 .
figure17.In Phase A (figures 17(a) and (b)), the saturation is unattainable during the simulation because the simulation collapses due to a numerical problem at large amplitude.We scan β h for this n = 1 case, whose results are shown in figure18.Basically, it shows that it is MHD unstable with γ ≈ 0.6% at β h = 0 (Note that the finite frequency when β h = 0% arises from the toroidal rotation).This may explain the absence of

Figure 15 .
Figure 15.Effect of β h on the frequency jump phenomenon showing time evolution of the mode frequency when β h = (a) 4.00%, (b) 3.14%, (c) 2.40%, and (d) 1.95%.A decrease in β h results in the delay of onset of the frequency jump or complete disappearance when β h is sufficiently low.

Figure 16 .
Figure 16.Toroidal rotation effect on the linear growth rate at two different experimental conditions studied in this paper.n represents the toroidal mode number and the growth rate is normalized to the Alfvén frequency ω A .

Figure 17 .
Figure 17.Evolution of the kinetic energy (E k ) and the mode frequency for the n = 1 mode at Phases A [(a), (b)] and C [(c), (d)], respectively.

Figure 18 .
Figure 18.Normalized growth rate of the linear mode versus energetic particle beta (β h ).The finite positive growth rate when β h = 0 indicates the n = 1 mode MHD unstable.

Figure 19 .
Figure 19.Evolution of the kinetic energy (E k ) and the mode frequency when β h = 4.0% [(a), (b)] and 2.3% [(c), (d)], respectively.Frequency jump disappears as β h is reduced which may reveal the coherence of EPs and the nonlinear mode.

Figure 20 .
Figure 20.Scan of the toroidal rotation speed at the magnetic axis (ω ϕ0 ) for the n = 1 in Phase C. The growth rate increases firstly and gradually decreases with the increase in ω ϕ0 .

Table 1 .
Comparisons between frequencies in linear stages of experiments and simulations in Phases A and C. The linear frequencies of experiments refer to the starting values of the frequency chirping.Mode number Experiments (kHz) Simulations (kHz) phase A n = 1