The EPEC code is employed to model the q95 windows for n = 2 and n = 1 resonant magnetic perturbation (RMP)-induced edge localized mode (ELM) suppression in typical KSTAR H-mode discharges. The plasma equilibria used in the study are derived by rescaling the experimental plasma equilibrium in KSTAR discharge #18594 measured at time t = 6450 ms. The n = 2 and n = 1 ELM suppression windows predicted by EPEC are comparatively narrow in q95 (i.e., Δq950.1) and are such that (assuming a relative pedestal pressure reduction of 15% is required to trigger ELM suppression), even in the middle of the windows, there is a threshold RMP coil current of about 1–3 kA/turn that must be exceeded before ELM suppression occurs. The n = 2 and n = 1 ELM suppression windows calculated by the EPEC code are consistent with the experimental observations and are also broadly similar to those recently calculated by the TM1 code [Hu et al., Phys. Plasmas 28, 052505 (2021)].

Magnetically diverted tokamak discharges operating in the high-confinement mode (H-mode)1 exhibit intermittent bursts of heat and particle transport, emanating from the outer regions of the plasma, which are known as (type-I) “edge localized modes” (ELMs).2 It is estimated that the heat load that ELMs will deliver to the tungsten plasma-facing components in a reactor-scale tokamak, such as ITER, will be large enough to cause a massive tungsten ion influx into the plasma core, and that the erosion associated with this process will unacceptably limit the lifetimes of these components.3 Consequently, the development of robust and effective methods for ELM control is a high priority for the international magnetic fusion program.

The most promising method for the control of ELMs in H-mode tokamak discharges is via the application of static “resonant magnetic perturbations” (RMPs). A complete RMP-induced ELM suppression was first demonstrated on the DIII-D tokamak.4 Subsequently, either mitigation or complete suppression of ELMs has been demonstrated on the JET,5 ASDEX-U,6 KSTAR,7 MAST,8 and EAST9 tokamaks.

RMP-induced ELM suppression is only observed to occur when q95 [the safety-factor value at the magnetic flux-surface that encloses 95% of the poloidal magnetic flux enclosed by the last closed magnetic flux-surface (LCFS)] lies in well-separated, narrow (i.e., Δq950.1) windows.10 Of course, other factors, such as plasma collisionality and toroidal rotation, also play a role in achieving the RMP-induced ELM suppression.11 However, in this paper, we shall concentrate on the role of q95. Recent simulations12–15 performed using the TM1 cylindrical, two-fluid, reduced-MHD (magnetohydrodynamics), initial value code,16–18 combined with the GPEC toroidal, ideal-MHD response code,19 have been able to account for virtually all salient features of the RMP-induced ELM suppression windows in q95 observed in the DIII-D and KSTAR tokamaks on the plausible assumption that the suppression is associated with a magnetic island chain driven at the top (i.e., the inner radius) of the pedestal region20 (i.e., the region of strong pressure and current gradients located at the edge of a typical H-mode tokamak discharge) that is sufficiently wide to reduce the plasma pressure just inside the inner radius of the pedestal by approximately 15%.

The TM1 code has two main deficiencies. First, it assumes a cylindrical geometry, with circular magnetic flux-surfaces, which is quite different from the toroidal, strongly shaped flux-surfaces typically found in H-mode tokamak plasmas. However, this deficiency is partly offset by the use of the fully toroidal GPEC code to determine the ideal response of the plasma to the applied RMP. Second, the TM1 code models RMP-induced changes in plasma rotation as being predominately poloidal in nature, when they are, in fact, predominately toroidal (because changes in the poloidal rotation are strongly suppressed by neoclassical poloidal flow damping). However, this deficiency is partly offset by artificially increasing the perpendicular momentum diffusivity by a factor of (Bϕ/Bθ)2 (which leads to the correct prediction for the mode locking threshold on the assumption that RMP-induced changes in the poloidal rotation are completely suppressed by the neoclassical flow damping21) Given the deficiencies in the TM1 code, it is desirable to obtain an independent confirmation of its predictions. This goal can be achieved by means of the EPEC code.22–24 

The EPEC code22–24 employs an asymptotic matching21,25–34 approach to determine the resistive response of a toroidal tokamak equilibrium to an applied RMP. According to this approach, the response of the plasma to the applied RMP is governed by a combination of flux-freezing and perturbed force balance—this combination is usually referred to as “marginally stable ideal-magnetohydrodynamics (MHD)”—everywhere in the plasma apart from a number of relatively narrow (in the radial direction) regions in which the applied perturbation resonates with the equilibrium magnetic field. Magnetic reconnection can take place within the resonant regions to produce relatively narrow magnetic island chains. Within the resonant regions, the plasma response is governed by linear or nonlinear two-fluid resistive-MHD, depending on whether the induced magnetic island widths are smaller or larger, respectively, than the corresponding linear layer widths. Thus, when employing the asymptotic matching approach, the equations of marginally stable ideal-MHD are solved in the so-called “outer region” that comprises most of the plasma (and the surrounding vacuum), the equations of linear/nonlinear two-fluid resistive-MHD are solved in the various resonant layers that constitute the so-called “inner region,” and the two sets of solutions are then asymptotically matched to one another. The EPEC model is described in detail in the  Appendix. (Note that, unlike the TM1 model, the EPEC model is fully toroidal and takes both poloidal and toroidal RMP-induced changes in plasma rotation into account.)

The main advantage of the asymptotic matching approach is that it removes the very short (and largely irrelevant) Alfvén time from the problem. In fact, the EPEC code is capable of accurately simulating the resistive response of a toroidal tokamak plasma to an RMP while taking time steps that extend over many Alfvén times. In this manner, the code is able to simulate a complete plasma discharge in a matter of minutes of real time.

The aim of this paper is to confirm the properties of the q95 windows for n = 2 (here, n is the toroidal mode number of the applied RMP) and n = 1 RMP-induced ELM suppression in typical KSTAR H-mode discharges that were previously calculated by the TM1 code, and are displayed in Fig. 8 of Ref. 15.

The Korea Superconducting Tokamak Advanced Research (KSTAR) device is a medium-sized tokamak that is equipped with a set of in-vessel magnetic field coils whose purpose is to generate an RMP capable of suppressing ELM crashes.35 There are 12 coils located on the outboard side of the device, and consisting of four equally spaced (in toroidal angle) top coils, four equally spaced middle coils, and four equally spaced bottom coils. Each coil has two turns. The field coils are capable of generating an RMP of toroidal mode number n = 1 or n = 2.

KSTAR discharge #18594 is a typical H-mode discharge in which an n = 2 RMP was employed to suppress ELMs.36,37

Figure 1 shows the (kinetic) plasma equilibrium in KSTAR discharge #18594 at time t = 6450 ms, at which time B0=1.79T,R0=1.80m,r100=0.598m, q95=4.04 (these quantities are defined in Appendix 2 a and 2 b), and ELM suppression has not yet occurred.

FIG. 1.

Left panel: contours of the equilibrium poloidal magnetic flux in KSTAR discharge #18594 at time t = 6450 ms. The scale major radius is R0=1.80m. The white dot indicates the magnetic axis, the white curve indicates the LCFS, and the thick black line indicates the limiter. Upper-right panel: safety-factor profile in KSTAR discharge #18594 at time t = 6450 ms. Lower-right panel: total plasma pressure profile in KSTAR discharge #18594 at time t = 6450 ms.

FIG. 1.

Left panel: contours of the equilibrium poloidal magnetic flux in KSTAR discharge #18594 at time t = 6450 ms. The scale major radius is R0=1.80m. The white dot indicates the magnetic axis, the white curve indicates the LCFS, and the thick black line indicates the limiter. Upper-right panel: safety-factor profile in KSTAR discharge #18594 at time t = 6450 ms. Lower-right panel: total plasma pressure profile in KSTAR discharge #18594 at time t = 6450 ms.

Close modal

Figure 2 shows the number density, temperature, impurity ion toroidal angular velocity, and perpendicular diffusivity profiles in KSTAR discharge #18594 at time t = 6450 ms. The majority and impurity ion number density profiles are determined from the measured electron number density profile on the assumption that the majority ions are deuterium, the impurity ions are carbon-VI, and that Zeff=2.0 throughout the plasma (see Appendix 3 a). (Incidentally, Zeff is not directly measured in KSTAR. The value 2.0 is a best guess based on the measured stored thermal energy.) Note that the perpendicular energy and momentum diffusivities (see Appendix 6 and 8) have both been given the plausible values of 1m2s1, whereas the perpendicular particle diffusivity has been given the plausible value 0.20m2s1. It can be seen that the top of the pedestal corresponds to ΨN=0.945. (The normalized poloidal flux ΨN is defined in Appendix 2 b.)

FIG. 2.

Top panel: the red, green, blue, cyan, and magenta curves show the electron number density (1019m3), electron temperature (keV), (thermal) ion temperature (keV), carbon-VI ion number density (1018m3), and flux-surface averaged neutral density (1017m3) (solid/dotted corresponds to high/low-recycling) profiles, respectively, in KSTAR discharge #18594 at time t = 6450 ms. Middle panel: carbon-VI impurity toroidal angular velocity profile (on outboard midplane) in KSTAR discharge #18594 at time t = 6450 ms. Bottom panel: the red, green, and blue curves show the perpendicular momentum (χϕ), energy (χE), and particle (D) diffusivity profiles, respectively, in KSTAR discharge #18594 at time t = 6450 ms. The common vertical dotted lines indicate the location of the top of the pedestal, ΨN=0.945.

FIG. 2.

Top panel: the red, green, blue, cyan, and magenta curves show the electron number density (1019m3), electron temperature (keV), (thermal) ion temperature (keV), carbon-VI ion number density (1018m3), and flux-surface averaged neutral density (1017m3) (solid/dotted corresponds to high/low-recycling) profiles, respectively, in KSTAR discharge #18594 at time t = 6450 ms. Middle panel: carbon-VI impurity toroidal angular velocity profile (on outboard midplane) in KSTAR discharge #18594 at time t = 6450 ms. Bottom panel: the red, green, and blue curves show the perpendicular momentum (χϕ), energy (χE), and particle (D) diffusivity profiles, respectively, in KSTAR discharge #18594 at time t = 6450 ms. The common vertical dotted lines indicate the location of the top of the pedestal, ΨN=0.945.

Close modal

Regrettably, there are no KSTAR data regarding neutrals inside the last closed magnetic flux-surface (LCFS). However, there are some DIII-D data regarding such neutrals.38 Given the similarity between KSTAR and DIII-D H-mode plasmas (in terms of size, shape, profiles, et cetera), it is not unreasonable to use the DIII-D data as a baseline with which to characterize KSTAR discharge #18594. Consequently, the flux-surface-averaged neutral number density profile is assumed to take the form

nn(r)=nn(r100)1+(rr100)2/ln2.
(1)

(The flux-surface label r is defined in Appendix 2 a. r100 is the flux-surface label of the LCFS. The flux-surface average operator is defined in Appendix 3 e.) We shall consider two cases: a “high-recycling” case in which nn(r100)=1.3×1017m3 (which is similar to the inferred DIII-D value), and a “low-recycling” case in which nn(r100)=2.5×1016m3. The parameter ln is given the value ln=0.013m.38 The flux-surface neutral poloidal asymmetry parameter (see Appendix 3 e),

yn=nnB2nnB2,
(2)

is given the value 1.5,38 which implies a greater concentration of neutrals at the X-point. Finally, the flux-surface averaged deuterium-atom/deuterium-ion charge exchange rate constant takes the value39 

σvicx=4×1014m3s1.
(3)

Unfortunately, there are no usable poloidal rotation data in KSTAR discharge #18594. Hence, it is necessary to deduce the E×B frequency profile from a combination of the measured density and temperature profiles, the measured impurity ion toroidal angular velocity profile, and neoclassical theory. According to the neoclassical theory,40 in the absence of any source of parallel momentum or heat flux (other than that due to charge exchange with neutrals), the E×B frequency profile is related to the impurity ion toroidal angular velocity profile as follows (see Appendix 3 f),

ωE(r)=ωϕI(r)ωθI(r)ωI(r).
(4)

Here, ωϕI,ωθI, and ωI are the impurity ion toroidal rotation, poloidal rotation, and diamagnetic frequencies, respectively. Figure 3 shows the various components of ωE calculated at the n = 2 resonant surfaces in the pedestal region of KSTAR discharge #18594 at time t = 6450 ms. It can be seen that all three terms on the right-hand side of the previous equation make non-negligible contributions to the net E×B frequency. In the low-recycling case, there is an Er well in the pedestal (where the ωE curve dips below zero). However, in the high-recycling case, this well is eliminated via the action of charge exchange with neutrals.

FIG. 3.

Components of the E×B frequency calculated at the various n = 2 resonant surfaces in the pedestal of KSTAR discharge #18594 at time t = 6450 ms. The black, red, green, and blue curves correspond to ωE, ωϕI,ωI, and ωθI, respectively. The resonant surfaces correspond to m = 6, 7, 8, 9, and 10, respectively, in order from the left to the right. The vertical dotted lines indicate the location of the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

FIG. 3.

Components of the E×B frequency calculated at the various n = 2 resonant surfaces in the pedestal of KSTAR discharge #18594 at time t = 6450 ms. The black, red, green, and blue curves correspond to ωE, ωϕI,ωI, and ωθI, respectively. The resonant surfaces correspond to m = 6, 7, 8, 9, and 10, respectively, in order from the left to the right. The vertical dotted lines indicate the location of the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

Close modal

Figure 4 shows the natural frequencies (in the absence of the RMP) calculated at the n = 2 resonant surfaces in the pedestal region of KSTAR discharge #18594 at time t = 6450 ms (see Appendix 6 for the definition of a natural frequency). It can be seen that the natural frequency predicted by linear theory [see Eq. (A84)] is strongly in the electron direction (i.e., positive) at the bottom of the pedestal, attains a maximum value in the middle of the pedestal, before passing through zero at the top of the pedestal, and then reversing sign in the plasma core. In the high-recycling case, the natural frequency predicted by nonlinear theory [see Eq. (A85)] exhibits a similar behavior, except that it does not peak in the middle of the pedestal. The nonlinear natural frequency is in the electron direction throughout most of the pedestal because of the influence of charge exchange with neutrals which, when the neutrals are concentrated at the X-point (as is assumed to be the case), tends to shift the frequency in this direction.23 In the low-recycling case, this shift is not as strong, and the nonlinear natural frequency is consequently in the ion direction throughout most of the pedestal.

FIG. 4.

Natural frequencies (in the absence of the RMP) calculated at the various n = 2 resonant surfaces in the pedestal of KSTAR discharge #18594 at time t = 6450 ms. The red, blue, and black curves show the natural frequencies predicted by linear theory, nonlinear theory, and the assumption that magnetic island chains are convected at the local E×B velocity, respectively. The resonant surfaces correspond to m = 6, 7, 8, 9, and 10, respectively, in order from the left to the right. The vertical dotted lines indicates the location of the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

FIG. 4.

Natural frequencies (in the absence of the RMP) calculated at the various n = 2 resonant surfaces in the pedestal of KSTAR discharge #18594 at time t = 6450 ms. The red, blue, and black curves show the natural frequencies predicted by linear theory, nonlinear theory, and the assumption that magnetic island chains are convected at the local E×B velocity, respectively. The resonant surfaces correspond to m = 6, 7, 8, 9, and 10, respectively, in order from the left to the right. The vertical dotted lines indicates the location of the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

Close modal

In order to search for q95 windows in which RMP-induced ELM suppression is possible, we require a collection of otherwise similar plasma equilibria with a range of different q95 values. Such a collection is obtained by rescaling the equilibrium discussed in Sec. II. The rescaling is performed in two steps.41 First, the original equilibrium (labeled “ori”) is rescaled to an intermediate equilibrium (labeled “int”). Second, the intermediate equilibrium is rescaled to a final equilibrium (labeled “fin”). The first rescaling is such that

gint(r)=[gori2(r)+α]1/2,
(5)

where

α=[q95target2qori2(r95)1]gori2(r95).
(6)

Here, R0B0g(r)ϕ is the toroidal magnetic field-strength, where ϕ is a geometric toroidal angle (see Appendix 2 b). Moreover, r95 is the value of the flux-surface label r at the surface at which ΨN=0.95 (see Appendix 2 b). The second rescaling is such that

Ψpfin(r)=βΨpori(r),
(7)
gfin(r)=βgint(r),
(8)
nafin(r)=βnaori(r),
(9)
Tafin(r)=βTaori(r),
(10)

where

β=gori(r100)gint(r100).
(11)

Here, the species label a stands for either e (electrons), i (majority ions), or I (impurity ions) (see Appendix 3 a). Moreover, Ψp(r) is the poloidal magnetic flux (divided by 2π) (see Appendix 2 b), na(r) is the species-a number density, and Ta(r) is the species-a temperature (in energy units). The net result of the rescaling procedure is a plasma equilibrium in which the toroidal plasma current has been modified in such a manner that q95 is changed from its original value of 4.04 to the new value q95target, while keeping the vacuum toroidal magnetic field-strength constant.

Obviously, the rescaling method used in the paper runs the risk of increasing the pedestal pressure gradient to such an extent that the peeling-ballooning stability threshold is violated.42 A more physical approach would be to incorporate the well-known EPED edge stability model43 into the rescaling method, in order to avoid this possibility. This will be attempted in future publications.

According to Eq. (A4), the ideal (i.e., Ψk=0 for all k) response of the plasma at the kth resonant surface to the applied RMP is governed by

ΔΨk=|Ekk|χk.
(12)

Here, ΔΨk is a (complex) measure of the strength of the ideal current-sheet that develops at the resonant surface, in response to the applied RMP. Moreover, Ekk is the (real) normalized tearing stability index at the resonant surface,25 and χk is a (complex) parameter determined by the GPEC code (see Appendix 2 d) that measures the ideal response of the kth resonant surface to the currents flowing in the RMP coils. In fact, the previous equation can also be written

ΔΨk=|Ekk|Irmp(einΔϕχkU+χkM+e+inΔϕχkL).
(13)

Here, Irmp is the (real) peak current flowing in the upper, middle, and lower coil sets (in units of kiloamps per turn). Moreover, χkU,χkM, and χkL are the (complex) χk values determined from GPEC when 1 kA/turn flows in each alone of the upper, middle, and lower coil sets, respectively, with the same toroidal phase. Finally, Δϕ is the (real) toroidal phase-shift between the current patterns flowing in the middle and upper, and also between the current patterns flowing in the lower and middle, coil sets, respectively. {In other words, if the nth toroidal harmonic of the current pattern flowing in the middle coil set varies with toroidal angle as cos(nϕ), then that in the upper coil set varies as cos[n(ϕΔϕ)], and that in the lower coil set as cos[n(ϕ+Δϕ)].} In this paper, the peak currents flowing in each coil set are assumed to be equal to one another, in accordance with standard experimental practice in KSTAR.44 Likewise, the toroidal phase-shift between the currents flowing in the middle and the upper coil sets is assumed to be the same as that between the currents flowing in the lower and the middle coil sets.

As will become clear in Sec. V, in this paper we associate RMP-induced ELM suppression with the formation of an RMP-driven island chain at the top of the pedestal (i.e., ΨN=0.945). Suppose that the closest resonant surface to the top of the pedestal is the kth resonant surface. In order to facilitate ELM suppression, we can optimize the value of the toroidal phase-shift between the current patterns flowing in the three coil sets so as to maximize the value of |ΔΨk|. (In fact, we maximize the modulus of a linear superposition of the ΔΨk values of the two resonant surface that straddle the top of the pedestal, with superposition weights that are inversely proportional to the distances of the resonant surfaces from the pedestal top.) If

χkUχkM=|χkUχkM|eiγkUM,
(14)
χkUχkL=|χkUχkL|eiγkUL,
(15)
χkMχkL=|χkMχkL|eiγkML,
(16)

where γkUM,γkUL, and γkML are real, then the optimization procedure boils down to maximizing the function

f(Δϕ)=|χkUχkM|cos(nΔϕγkUM)+|χkUχkL|cos(2nΔϕγkUL)+|χkMχkL|cos(nΔϕγkML).
(17)

This goal is achieved in the EPEC code numerically.

Figure 5 shows the optimum toroidal phase-shift, as determined by EPEC, for both n = 2 and n = 1 RMPs as a function of q95. (The optimal shift is independent of the recycling level.) It can be seen that, on average, the optimal phase-shift for n = 2 RMPs is 90°, which is in accordance with KSTAR experimental observations.37 However, the actual optimal phase-shift is slightly less than 90° for large q95 (i.e., q95>3.7), and slightly larger for small q95. The optimal phase-shift for n = 1 RMPs exhibits a much greater variation with q95. The phase-shift used in experiments is 270°.44 It can be seen that the optimal phase-shift determined by EPEC approaches this value at small q95, but is significantly different at large q95.

FIG. 5.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 or an n = 1 RMP. Top panel: optimum toroidal phase-shift between the middle and top, and also between the bottom and middle, RMP coil sets for an n = 2 RMP. Bottom panel: optimum toroidal phase-shift between the middle and top, and also between the bottom and middle, RMP coil sets for an n = 1 RMP.

FIG. 5.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 or an n = 1 RMP. Top panel: optimum toroidal phase-shift between the middle and top, and also between the bottom and middle, RMP coil sets for an n = 2 RMP. Bottom panel: optimum toroidal phase-shift between the middle and top, and also between the bottom and middle, RMP coil sets for an n = 1 RMP.

Close modal

Figures 6 and 8 show the results obtained from an EPEC simulation of the response of a typical KSTAR H-mode plasma to an n = 2 RMP. (The compete EPEC model is outlined in the  Appendix.) An n = 2 current of peak amplitude Irmp=4kA/turn is applied to each row of RMP coils with the optimal phasing. Using the rescaled plasma equilibria described in Sec. III, q95 is steadily ramped downward from a value of about 5 to a value of about 3 over a 2.2 s time interval. Note that all resonant surfaces in the region 0.995<ΨN<1.00 are neglected in the calculation (because GPEC does not give reliable results for such surfaces).

FIG. 6.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=1.3×1017m3. Top panel: n = 2 natural frequencies, in the absence of RMP, as functions of q95. Bottom panel: n = 2 natural frequencies, in the presence of RMP, as functions of q95. The black, red, green, blue, yellow, cyan, magenta, brown, and pink curves correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively.

FIG. 6.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=1.3×1017m3. Top panel: n = 2 natural frequencies, in the absence of RMP, as functions of q95. Bottom panel: n = 2 natural frequencies, in the presence of RMP, as functions of q95. The black, red, green, blue, yellow, cyan, magenta, brown, and pink curves correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively.

Close modal

Figures 6 and 7 show the natural frequencies, both in the absence and in the presence of the applied RMP, of the various n = 2 resonant surfaces present in the plasma as functions of q95, in the high- and low-recycling cases, respectively. In both cases, as a resonant surface passes through the pedestal (from the inner to the outer radius of the pedestal—recall that q95 is ramped down in the simulations), its natural frequency starts off being comparatively small and in the ion direction (i.e., negative) in the plasma core, passes through zero close to the top of the pedestal, peaks in the center of the pedestal, and ends up being comparatively large and in the electron direction (i.e., positive) at the edge of the pedestal. It can be seen that if the magnitude of the natural frequency in the absence of the RMP falls below about 5krad/s, then the mode locks to the RMP (i.e., its natural frequency in the presence of the RMP becomes zero). Note that the differences between the low- and high-recycling cases are comparatively minor.

FIG. 7.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=2.5×1016m3. Top panel: n = 2 natural frequencies, in the absence of RMP, as functions of q95. Bottom panel: n = 2 natural frequencies, in the presence of RMP, as functions of q95. The black, red, green, blue, yellow, cyan, magenta, brown, and pink curves correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively.

FIG. 7.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=2.5×1016m3. Top panel: n = 2 natural frequencies, in the absence of RMP, as functions of q95. Bottom panel: n = 2 natural frequencies, in the presence of RMP, as functions of q95. The black, red, green, blue, yellow, cyan, magenta, brown, and pink curves correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively.

Close modal

Figures 8 and 9 show the full widths (in ΨN) of the magnetic island chains induced in the plasma by the applied RMP, both in the absence and in the presence of shielding due to plasma flow, as functions of q95, in the high- and low-recycling cases, respectively. In the absence of shielding, the so-called “vacuum” island chains are all wide, extend over most of the pedestal, and would cause a very significant reduction in the pedestal pressure. In the presence of shielding, the island chains are all narrow and have a negligible effect on the pedestal pressure. The only exceptions to this state of affairs occur when one of the island chains locks to the RMP. This happens when the chain's resonant surface passes through the top of the pedestal (because its associated natural frequency simultaneously passes through zero) and leads to the formation of a relatively wide (ΔΨN0.1) magnetic island chain at the top of the pedestal. Again, the differences between the low- and high-recycling cases is comparatively minor.

FIG. 8.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=1.3×1017m3. Top panel: full n = 2 vacuum island widths as functions of q95. Bottom panel: full n = 2 island widths as functions of q95. The black, red, green, yellow, cyan, magenta, brown, and pink areas correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively. The horizontal dotted lines indicate the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995.

FIG. 8.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=1.3×1017m3. Top panel: full n = 2 vacuum island widths as functions of q95. Bottom panel: full n = 2 island widths as functions of q95. The black, red, green, yellow, cyan, magenta, brown, and pink areas correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively. The horizontal dotted lines indicate the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995.

Close modal
FIG. 9.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=2.5×1016m3. Top panel: full n = 2 vacuum island widths as functions of q95. Bottom panel: full n = 2 island widths as functions of q95. The black, red, green, yellow, cyan, magenta, brown, and pink areas correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively. The horizontal dotted lines indicate the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995.

FIG. 9.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP, assuming that nn(r100)=2.5×1016m3. Top panel: full n = 2 vacuum island widths as functions of q95. Bottom panel: full n = 2 island widths as functions of q95. The black, red, green, yellow, cyan, magenta, brown, and pink areas correspond to m = 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively. The horizontal dotted lines indicate the top of the pedestal, ΨN=0.945, as well as the effective edge of the plasma, ΨN=0.995.

Close modal

Note that the EPEC model fully takes into account the strong poloidal coupling of tearing modes with different poloidal mode numbers, but a common toroidal mode number, due to the Shafranov shifts and departures from circular cross sections of magnetic flux-surfaces, that is predicted by linear MHD theory.30 This coupling manifests itself in the model via the off diagonal elements of the Ekk matrix (see Appendix 2 d). However, it turns out that the differential flow in the pedestal of KSTAR H-mode plasmas is so strong that substantial driven reconnection is only possible at one particular resonant surface in the plasma (i.e., a surface whose natural frequency is comparatively close to zero). In other words, the differential flow overwhelms the MHD coupling and effectively decouples the various tearing modes. Consequently, the response of the plasma to the applied RMP is similar, in many respects, to that predicted by cylindrical theory (which completely neglects the MHD coupling). This fact may account for the success of the cylindrical TM1 code in predicting the response of H-mode plasmas to RMPs. The decoupling of tearing modes via differential flow accounts for the fact that the appearance of a wide locked magnetic island chain at a resonant surface whose natural frequency is close to zero does not lead to substantial driven reconnection at neighboring surfaces, despite the strong MHD coupling between these surfaces.

Finally, Fig. 10 shows the cumulative change in the pedestal pressure due to the magnetic island chains induced in the pedestal by the applied RMP (see Appendix 8) as a fraction of the pressure at the top of the pedestal (i.e., ΨN=0.945), calculated as functions of q95 and the peak RMP coil current, Irmp. This figure was generated by combining 41 simulations of the form shown in Figs. 6–9, each performed with a different value of Irmp. It can be seen that the relative pressure reduction is negligible except when q95 lies in certain narrow (i.e., Δq950.1) windows that are evenly spaced in q95 with an interval of 0.5. The first, second, third, and fourth windows from the left correspond to comparatively wide m = 6, 7, 8, and 9 magnetic island chains, respectively, induced at the top of the pedestal. Note that the windows in the low-recycling case are slightly wider than those in the high-recycling case and also shifted downward, by about 0.05, in q95.

FIG. 10.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP. The plot shows the cumulative pressure change due to the n = 2 magnetic island chains induced in the plasma by the RMP as a fraction of the pressure at the top of the pedestal (i.e., ΨN=0.945), plotted as a function of q95 and the peak current flowing in the RMP coils, Irmp. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

FIG. 10.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 2 RMP. The plot shows the cumulative pressure change due to the n = 2 magnetic island chains induced in the plasma by the RMP as a fraction of the pressure at the top of the pedestal (i.e., ΨN=0.945), plotted as a function of q95 and the peak current flowing in the RMP coils, Irmp. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

Close modal

Now, the experimental m=8/n=2 RMP-induced ELM suppression window in KSTAR occurs at q953.8,45 rather than 4.0 (as indicated by the figure); moreover, the experimental m=7/n=2 ELM suppression window occurs at 3.4 rather than 3.5.44 However, q95 in KSTAR experiments is determined from a magnetic plasma equilibrium (in which the influence of the bootstrap current and pressure on the equilibrium is neglected), whereas q95 in EPEC is determined from a kinetic plasma equilibrium (in which the influence of the bootstrap current and pressure is taken into account). It turns out that q95 for a magnetic equilibrium is smaller than that for a corresponding kinetic equilibrium by about 0.2. Hence, there is, in fact, reasonably a good agreement between the EPEC predictions of locations of the n = 2 q95 windows and experimental observations. (It should be noted that the q95 windows in Fig. 8 of Ref. 15 have been shifted down by 0.2 in order to take into account the difference between the q95 values obtained from magnetic and kinetic equilibria.) The main uncertainty in the EPEC model is associated with the poorly known neutral particle profile within the LCFS. As we have seen, charge exchange with neutrals can profoundly affect the E×B velocity profile in the pedestal, which, in turn, leads to a shift in the locations of the q95 windows. Unfortunately, data on the neutral particles within the LCFS of H-mode tokamak plasmas are presently very hard to come by.

The principal factors that affect the widths (in q95) of the ELM suppression windows are the kinetic equilibrium (which determines the χk factors calculated by GPEC), the shear in the natural frequency close to the top of the pedestal (which depends on the Er, number density, and electron/ion temperature profiles in this region), the plasma resistivity close to the top of the pedestal (which depends on the number density, electron temperature, and Zeff values in this region), and the neoclassical poloidal viscosity close to the top of the pedestal (which depends on the number density, ion temperature, and Zeff values in this region). The widths are far less sensitive to factors, such as the momentum diffusivity and the tearing stability matrix, Ekk. Although it is hard to put a number on the uncertainty in the widths, we note that all of the quantities upon which the widths principally depend ought to be relatively well-determined in the discharge under consideration.

If we accept the hypothesis of Ref. 15 that a relative reduction in the pedestal pressure of 15% is required to trigger ELM suppression, then it is clear that there is a threshold RMP coil current that must be exceeded in order to trigger ELM suppression, even when q95 lies in the middle of one of the ELM suppression windows. This threshold value is about 1–3 kA/turn and increases with the increasing q95. Note that the threshold does not represent a threshold for locked island formation, but rather a threshold for the formation of a locked island whose width is sufficient to significantly affect the pedestal pressure profile. Hence, there is no contradiction between the results of this paper and those reported in Refs. 22–24. The threshold was not manifest in the latter papers because they did not calculate the reduction in the pedestal pressure due to the island chain.

If we compare the n = 2 ELM suppression windows shown in Fig. 10 with the corresponding windows previously calculated by the TM1 code (and shown in Fig. 8 of Ref. 15), taking into account that the TM1 q95 windows have been shifted down by 0.2, then we can see that there is good agreement between the predictions of the EPEC and TM1 codes. In particular, the positions (in q95), widths (in q95), and heights (in Irmp) of the windows predicted by the two codes are about the same. (The TM1 simulations exhibit a monotonic decrease in the pressure reduction, from window to window, as q95 increases. The EPEC simulations exhibit a similar trend, but in a more nuanced fashion, in that the decrease is not quite monotonic. This slight difference is probably due to the additional neoclassical and neutral density physics incorporated into EPEC.) Incidentally, as is clear from Fig. 8 of Ref. 15, there is a good agreement between the n = 2 ELM suppression windows predicted by the TM1 code and those observed experimentally in KSTAR.

Figure 11 shows a similar calculation to Fig. 10, except that an n = 1 RMP is applied to the plasma. The RMP is generated by running n = 1 currents through the three rows of KSTAR RMP coils with the optimal phasing. As before, Fig. 11 is produced by combining the simulations similar in form to those shown in Figs. 6–9, each performed with a different value of Irmp. In each simulation, q95 is ramped downward from a value of about 5.5 to a value of about 2.5 over a time interval of 3 s.

FIG. 11.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 1 RMP. The plot shows the cumulative pressure change due to the n = 1 magnetic island chains induced in the plasma by the RMP as a fraction of the pressure at the top of the pedestal (i.e., ΨN=0.945), plotted as a function of q95 and the peak current flowing in the RMP coils, Irmp. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

FIG. 11.

EPEC simulation of the response of set of self-similar discharges with a range of different q95 values that are derived from KSTAR discharge #18594 at time t = 6450 ms to an n = 1 RMP. The plot shows the cumulative pressure change due to the n = 1 magnetic island chains induced in the plasma by the RMP as a fraction of the pressure at the top of the pedestal (i.e., ΨN=0.945), plotted as a function of q95 and the peak current flowing in the RMP coils, Irmp. The top panel corresponds to nn(r100)=1.3×1017m3, whereas the bottom panel corresponds to nn(r100)=2.5×1016m3.

Close modal

It can be seen that the relative reduction in the pedestal pressure due to the island chains induced by the RMP is negligible except when q95 lies in certain narrow (i.e., Δq950.1) windows that are evenly spaced in q95 with an interval of 1.0. The first, second, and third windows from the left correspond to comparatively wide m = 3, 4, and 5 magnetic island chains, respectively, induced at the top of the pedestal. As before, the windows in the low-recycling case are slightly wider than those in the high-recycling case and also shifted downward, by about 0.05, in q95.

If we again accept the hypothesis of Ref. 15 that a relative reduction in the pedestal pressure of 15% is required to trigger ELM suppression, then it is clear that there is a threshold RMP coil current that must be exceeded in order to trigger ELM suppression, even when q95 lies in the middle of one of the ELM suppression windows. This threshold value is about 1–3 kA/turn and increases with the increasing q95.

If we compare the n = 1 ELM suppression windows shown in Fig. 11 with the corresponding windows previously calculated by the TM1 code (and shown in Fig. 8 of Ref. 15), taking into account that the TM1 q95 windows have been shifted down by 0.2, then we can see that there is a good agreement between the predictions of the EPEC and TM1 codes. (As before, as is clear from Fig. 8 of Ref. 15, there is good agreement between the n = 1 ELM suppression windows predicted by the TM1 code and those observed experimentally in KSTAR.) The widths (in q95) and heights (in Irmp) of the windows predicted by the TM1 code are somewhat larger than those predicted by the EPEC code. On the other hand, the locations of the windows (in q95) are about the same. As before, the main uncertainty in the EPEC model is associated with the poorly known neutral particle profile within the LCFS.

It should be noted that the application of n = 1 RMPs to H-mode tokamak plasmas can be hazardous, because the RMPs can trigger disruption-producing locked modes (in particular, those associated with the m=2/n=1 mode). Hence, it is advisable to optimize the RMP spectrum in such a manner that the RMP resonates strongly with tearing modes whose resonant surfaces lie close to the top of the pedestal without simultaneously resonating strongly with the tearing modes whose resonant surfaces lie in the plasma core.46 

We have employed the EPEC code (see the  Appendix) to model the q95 windows for n = 2 and n = 1 RMP-induced ELM suppression in typical KSTAR H-mode discharges. The plasma equilibria used in this study are derived by rescaling (see Sec. III) the experimental (kinetic) plasma equilibrium in KSTAR discharge #18594 measured at time t = 6450 ms (see Sec. II).

Our findings are summarized in Figs. 10 and 11, which show the n = 2 and n = 1 ELM suppression windows, respectively, in q95 vs Irmp space. Here, q95 is the safety-factor at the 95% magnetic flux-surface, whereas Irmp is the peak current applied to the three rows of RMP coils in KSTAR. It can be seen that the ELM suppression windows are comparatively narrow in q95 (i.e., Δq950.1) and that (assuming a relative pedestal pressure reduction of 15% is required to trigger ELM suppression15) even in the middle of the windows, there is a threshold RMP coil current of about 1–3 kA/turn that must be exceeded before ELM suppression occurs.

The KSTAR n = 2 and n = 1 ELM suppression windows calculated by the EPEC code are broadly similar to those previously calculated by the TM1 code, and shown in Fig. 8 of Ref. 15. (Note that the TM1 windows have been shifted downward in q95 by 0.2 in order to take into account the difference between the magnetic equilibria used experimentally and the kinetic equilibria used in the TM1 code.) The small differences [in terms of width (in q95), height (in Irmp), and position (in q95)] between the ELM suppression windows calculated by the EPEC and TM1 codes is probably related to the different (and largely complementary) sets of approximations made in the two codes. For instance, the EPEC code is fully toroidal, whereas the TM1 code is cylindrical; the EPEC code interpolates between the linear and the nonlinear resonant response regimes, whereas the TM1 code calculates the resonant response directly; the EPEC code allows for warm ions, whereas the TM1 code employs a cold-ion approximation; the EPEC code allows for both poloidal and toroidal RMP-induced plasma angular velocity changes, whereas the TM1 code represents toroidal angular velocity changes as poloidal angular velocity changes (with a rescaled perpendicular momentum diffusivity to get the correct mode locking threshold); the EPEC code uses analytic theory to determine the temperature and density flattening across each magnetic island chain induced in the plasma by the RMP (taking into account that the parallel transport is convective, rather than diffusive, in nature), whereas the TM1 code solves transport equations self-consistently to determine the flattening (but models the parallel transport as diffusive).

The main uncertainty in the EPEC model is associated with the poorly known neutral particle profile within the LCFS. Charge exchange with neutrals can profoundly affect the E×B velocity profile in the pedestal, which, in turn, leads to a shift in the locations of the q95 windows. Unfortunately, data on neutral particles within the LCFS of H-mode tokamak plasmas are presently very hard to come by.

The authors would like to thank Y. In for facilitating the research presented in this paper. The authors would also like to thank J.-K. Park and N. C. Logan for their advice on how to run the GPEC code. Part of the data analysis was performed using the OMFIT integrated modeling framework.47 

This research was supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, under Contract Nos. DE-FG02–04ER54742, DE-SC0021156, and DE-SC0020372. This research was also supported by R&D Program of “KSTAR Experimental Collaboration and Fusion Plasma Research (EN2101–12)” through the Korea Institute of Fusion Energy (KFE) funded by government funds.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1. Introduction

The Extended Perturbed Equilibrium Code (EPEC) model was introduced in Ref. 22, but slightly modified and corrected in Refs. 23 and 24. To avoid confusion, the latest (and, hopefully, final) version of the complete model is summarized in this  Appendix.

The main changes made to the model are as follows. In Eqs. (A8), (A81), and (A83), more accurate (i.e., fully toroidal) expressions have been substituted for those employed previously. In Eqs. (A53) and (A86), a charge exchange term that was erroneously omitted in previous versions of the model has been added. The neoclassical theory of impurity ion rotation has been incorporated into the model in Appendix 3 f. The determination of the linear layer width has been greatly improved (see Apppendix 5). Charge exchange damping has been incorporated into the plasma angular velocity evolution equations (see Appendix 6). The crucial model expression (A87), for the natural frequency as a function of the island width has been changed somewhat from that employed in previous studies. Finally, the model for the effect of induced magnetic islands on the plasma pressure profile has been greatly improved (see Appendix 8).

2. Plasma response in outer region
a. Coordinates

Let R, ϕ, and Z be right-handed cylindrical coordinates whose symmetry axis corresponds to the toroidal symmetry axis of the plasma. Let r, θ, and ϕ be right-handed flux coordinates whose Jacobian is J(r×θ·ϕ)1=rR2/R0. Here, R0 is a convenient scale major radius, r is a magnetic flux-surface label with dimensions of length, and θ is an axisymmetric angular coordinate that increases by 2π radians for every poloidal circuit of the magnetic axis. Let r = 0 correspond to the magnetic axis, and let r=r100 correspond to the LCFS. Let θ = 0 correspond to the inboard midplane, and let 0<θ<π correspond to the region above the midplane.

b. Equilibrium magnetic field

The equilibrium magnetic field is written B=R0B0[f(r)ϕ×r+g(r)ϕ], where B0 is a convenient scale toroidal magnetic field-strength, and q(r)=rg/(R0f) is the safety-factor profile. The equilibrium poloidal magnetic flux (divided by 2π), Ψp(r), satisfies dΨp/dr=R0B0f(r), where, by convention, Ψp(r100)=0. The normalized poloidal magnetic flux, ΨN(r), is defined such that ΨN(r)=1Ψp(r)/Ψp(0). Hence, ΨN(0)=0 and ΨN(r100)=1. Finally, if ΨN(r95)=0.95, then q95q(r95).

c. Perturbed magnetic field

Consider the response of the plasma to an RMP with n > 0 periods in the toroidal direction. We can write the components of the perturbed magnetic field in the form22 

rR2δB·rR02=ijψj(r)ei(mjθnϕ),
(A1)
JδB·ϕ×r=jΞj(r)ei(mjθnϕ),
(A2)
R2δB·ϕ=njΞj(r)mjei(mjθnϕ),
(A3)

where the sum is over all relevant poloidal harmonics of the perturbed magnetic field.

Let there be K resonant magnetic flux-surfaces in the plasma, labeled 1 through K. Consider the kth resonant surface, r = rk, at which nq(rk)=mk, where mk is a positive integer. Let Ψk=ψk(rk)/mk, and ΔΨk=[Ξk]rkrk+.22 Here, Ψk is the (complex) reconnected helical magnetic flux (divided by 2πR0) at the kth resonant surface, whereas ΔΨk (which has the same dimensions as Ψk) is a (complex) measure of the strength of the current sheet (consisting of filaments running parallel to the local equilibrium magnetic field) at the same resonant surface.

d. Toroidal tearing mode dispersion relation

In the presence of the RMP, the Ψk and the ΔΨk values are related according to the inhomogeneous toroidal tearing mode dispersion relation, which takes the form21,22,29–34

ΔΨk=k=1,KEkkΨk+|Ekk|χk.
(A4)

Here, Ekk (for k, k=1, K) is the dimensionless, Hermitian, toroidal tearing mode stability matrix, whereas the χk (for k=1,K) parameterize the current sheets driven at the various resonant surfaces when the plasma responds to the applied RMP in accordance with the equations of linearized, marginally stable, ideal-MHD.

The EPEC model determines the elements of the Ekk matrix using a high-q approximation. In fact, if Fkk is the inverse of the Ekk matrix, then22 

Fkk=G(Rk,Zk;Rk,Zk)ei(mkθkmkθk)dθk2πdθk2π
(A5)

and

G(Rk,Zk;Rk,Zk)=(1)nπ2RkRk/R02Γ(1/2)Γ(n+1/2)[cosh ηkkRk2+Rk2+(ZkZk)2]1/2×[(n1/2)P1/2n1(cosh ηkk)+P1/2n+1(cosh ηkk)n+1/2],
(A6)

with

ηkk=tanh1[2RkRkRk2+Rk2+(ZkZk)2].
(A7)

Here, the double integral in Eq. (A5) is taken around the kth resonant surface (cylindrical coordinates Rk, 0, Zk; flux coordinates rk, θk, 0, with rk constant; resonant poloidal mode number mk) and the k th resonant surface (cylindrical coordinates Rk, 0, Zk; flux coordinates rk,θk, 0, with rk constant; resonant poloidal mode number mk). Finally, the Γ(z) and Pμν(z) are the gamma functions and associated Legendre functions, respectively.

The (complex) χk parameters are determined from the GPEC code.19,48 To be more exact, the GPEC code calculates the (complex) dimensionless Δmkn parameters, which are defined in Eq. (44) of Ref. 49, and which measure the strengths of the ideal current sheets that develop at the various resonant magnetic flux-surfaces in the plasma in response to the applied RMP. The Δmkn parameters are related to the χk parameters according to22 

χkR0B0=iΔmkn|Ekk|(rkR0)2g(rk)mk[akk(rk)+(rk/R0qk)2],
(A8)

where qk=mk/n, and akk(r)=|r|2dθ/(2π). Note that the previous equation is a slightly more accurate expression than Eq. (A76) of Ref. 23.

3. Neoclassical physics
a. Plasma species

The plasma is assumed to consist of three (charged) species, namely, electrons (e), majority ions (i), and impurity ions (I). The charges of the three species are ee=e, ei = e, and eI=ZIe, respectively, where e is the magnitude of the electron charge. Quasi-neutrality demands that ne=ni+ZInI, where na(r) is the species-a number density. Let αI(r)=ZI(Zeff1)/(ZIZeff), where Zeff(r)=(ni+ZI2nI)/ne is the effective ion charge number. It follows that ni/ne=(ZIZeff)/(ZI1) and nI/ne=(Zeff1)/[ZI(ZI1)].

b. Collisionality parameters

Consider an equilibrium magnetic flux-surface whose label is r. Let

1γ(r)=qgBR2B0R02dθ2π,
(A9)

where B=|B|. It is helpful to define a new poloidal angle Θ such that

dΘdθ=γqgBR2B0R02.
(A10)

Let

I1=B0BdΘ2π,
(A11)
I2=BB0dΘ2π,
(A12)
I3=(BΘ)21B0BdΘ2π,
(A13)
I4,j=2jcos(jΘ)B/B0dΘ2π,
(A14)
I5,j=2jcos(jΘ)2(B/B0)2dΘ2π,
(A15)
I6(λ)=1λB/BmaxB/B0dΘ2π,
(A16)

where Bmax is the maximum value of B on the magnetic flux-surface, and j a positive integer. The species-a transit frequency is written ωta(r)=KtγvTa, where

Kt(r)=I12I3I22j=1,I4,jI5,j,
(A17)

and vTa=2Ta/ma.40 Here, ma is the species-a mass, and Ta(r) the species-a temperature (in energy units). The fraction of circulating particles is40 

fc(r)=3I24B02Bmax201λdλI6(λ).
(A18)

Finally, the dimensionless species-a collisionality parameter is written40 νa(r)=Kgt/(ωtaτaa), where gt(r)=fc/(1fc),

K(r)=38πI2I3Kt2,
(A19)
1τaa(r)=43π4πnaea4lnΛ(4πϵ0)2ma2vTa3.
(A20)

Here, the Coulomb logarithm, lnΛ, is assumed to take the same large constant value (i.e., lnΛ17), independent of the species.

c. Collisional friction matrices

Let xab=vTb/vTa. In the following, all quantities that are of order (me/mi)1/2,(me/mI)1/2, or smaller, are neglected with respect to unity. The 2 × 2 dimensionless ion collisional friction matrices, [Fii](r),[FiI](r),[FIi](r), and [FII](r), are defined to have the following elements:40,50,51

F00ii=αI(1+mi/mI)(1+xiI2)3/2,
(A21)
F01ii=32αI(1+mi/mI)(1+xiI2)5/2,
(A22)
F11ii=2+αI[13/4+4xiI2+(15/2)xiI4](1+xiI2)5/2,
(A23)
F01iI=32TiTIαI(1+mI/mi)xiI(1+xIi2)5/2,
(A24)
F11iI=274TiTIαIxiI2(1+xiI2)5/2,
(A25)
F11Ii=274αIxiI2(1+xiI2)5/2,
(A26)
F11II=TiTI{2αI2xIi+αI[15/2+4xiI2+(13/4)xiI4](1+xiI2)5/2},
(A27)
F10ii=F01ii,F00iI=F00ii,F10iI=F01ii,F00Ii=F00ii,F01Ii=F01ii,F10Ii=F01iI,F00II=F00ii,F01II=F01iI,F10II=F01iI.

The 2 × 2 dimensionless electron collisional friction matrix, [Fee](r), is defined to have the following elements:40,50,51F00ee=Zeff, F01ee=(3/2)Zeff, F10ee=F01ee, F11ee=2+(13/4)Zeff.

d. Neoclassical viscosity matrices

The 2 × 2 dimensionless species-a neoclassical viscosity matrix, [μa](r), is defined to have the following elements: μ00a=K00a,μ01a=(5/2)K00aK01a, μ10a=μ01a, μ11a=K11a5K01a+(25/4)K00a.40,50,51 Here,

Kjke=gt43π0exx4+j+kνDe(x)dx[x2+νeνDe(x)][x2+(5π/8)(ωteτee)1νTe(x)],
(A28)
νDe=3π4[(112x)ψ(x)+ψ(x)]+3π4Zeff,
(A29)
νϵe=3π2[ψ(x)ψ(x)],
(A30)
νTa(x)=3νDa(x)+νϵa(x),
(A31)

and

ψ(x)=2π0xet2dt2πxex2,
(A32)
ψ(x)=2πxex2.
(A33)

Furthermore,

Kjki=gt43π0exx2+j+kνDi(x)dx[x+νiνDi(x)][x+(5π/8)(ωtiτii)1νTi(x)],
(A34)
νDi=3π4[(112x)ψ(x)+ψ(x)]1x+3π4αI[(1xiI2x)ψ(xxiI)+ψ(xxiI)]1x,
(A35)
νϵi=3π2[ψ(x)ψ(x)]1x+3π2αI[mimIψ(xxiI)ψ(xxiI)]1x,
(A36)

and, finally,

KjkI=gt43π0exx2+j+kνDI(x)dx[x+νIνDI(x)][x+(5π/8)(ωtIτII)1νTI(x)],
(A37)
νDI=3π4[(112x)ψ(x)+ψ(x)]1x+3π41αI[(1xIi2x)ψ(xxIi)+ψ(xxIi)]1x,
(A38)
νϵI=3π2[ψ(x)ψ(x)]1x+3π21αI[mImiψ(xxIi)ψ(xxIi)]1x.
(A39)

Note that our expressions for the neoclassical viscosity matrices interpolate in the most accurate manner possible between the three standard neoclassical collisionality regimes (i.e., the banana, plateau, and Pfirsch–Schüter regimes).40 

e. Parallel force and heat balance

Let [μ̃I]=αI2(Ti/TI)xIi[μI]. The requirement of equilibrium force and heat balance parallel to the magnetic field38,40,50,51 leads us to define four 2 × 2 dimensionless ion matrices, [Lii](r),[LiI](r),[LIi](r), and [LII](r), where

([Lii],[LiI][LIi][LII])=([Fii+μi+Yin/yn],[FiI][FIi],[FII+μ̃I])1([Fii+Yin],[FiI][FIi],[FII]),
(A40)

the additional four 2 × 2 dimensionless ion matrices, [Gii](r),[GiI](r),[GIi](r), and [GII](r), where

([Gii],[GiI][GIi][GII])=τiiσvicxnn([Fii+μi+Yin/yn],[FiI][FIi],[FII+μ̃I])1,
(A41)

and the 2 × 2 dimensionless electron matrix, [Qee](r), where [Qee]=[Fee+μe]1. Here,38 

[Yin]=τiiσvicxnn[1,00,En/Ti],
(A42)
yn=nnB2nnB2,
(A43)

where

A(r)A(r,Θ)dΘB(r,Θ)/dΘB(r,Θ).
(A44)

Moreover, σvicx is the flux-surface averaged rate constant for charge exchange reactions between neutrals and majority ions, nn(r,Θ) the neutral particle number density, and En/Ti the ratio of the incoming neutral energy to the majority ion energy. The parameter yn takes into account the fact that the incoming neutrals at the edge of an H-mode tokamak plasma are usually concentrated at the X-point (i.e., yn>1).38 

f. Impurity ion angular rotation velocities

Let

ωθI(r)=VI·θB·θR0B0gR2,
(A45)
ωϕI(r)=VI·ϕ,
(A46)

where VI is the impurity ion fluid velocity, and the right-hand sides are evaluated on the outboard midplane. According to neoclassical theory40 

ωθI=KθωncI,
(A47)
ωϕI=ωE+ωI+ωθI,
(A48)

where

Kθ(r)=R02B02g2R2B2,
(A49)
ωE(r)=dΘdΨp,
(A50)
ωa(r)=TaeadlnpadΨp,
(A51)
ηa(r)=dlnTadlnna,
(A52)

and

ωncI(r)=G00IiωE[L00IIL01II(ηI1+ηI)]ωI[L00IiL01Ii(ηi1+ηi)]ωi.
(A53)

Here, pa(r)=naTa, and Φ(r) is the equilibrium electric scalar potential. Note the presence of the first term on the right-hand side of the previous equation, which was erroneously omitted in Ref. 22. This term relates to the influence of charge exchange with neutrals on the neoclassical impurity ion poloidal rotation frequency.

4. Plasma response in inner region

Let Ψk=R0B0Ψ̂keiφk,χk=R0B0χ̂keiζk, and Ekk=Êkkeiξkk, where Ψ̂k>0,φk,χ̂k>0, ζk, Êkk>0, and ξkk are all real quantities. Furthermore, let Xk=Ψ̂kcosφk and Yk=Ψ̂ksinφk. The resonant plasma response model at the kth resonant surface takes the form24 

(Ŵk+δ̂k)Sk(dXkdt̂+ϖ̂kYk)=k=1,KÊkk(cosξkkXksinξkkYk)+Êkkχ̂kcosζk,
(A54)
(Ŵk+δ̂k)Sk(dYkdt̂ϖ̂kXk)=k=1,KÊkk(cosξkkYk+sinξkkXk)+Êkkχ̂ksinζk,
(A55)

where

Sk=τR(rk)τA,
(A56)
τA=[μ0ρ(0)r1002B02]1/2,
(A57)
τR(r)=μ0r2σeeQ00ee,
(A58)
σee(r)=nee2τeeme,
(A59)
Ŵk=2Iϵ100r̂k(qgs)rk1/2(Xk2+Yk2)1/4,
(A60)
s(r)=dlnqdlnr,
(A61)
δ̂k=δlinear(rk)R0ϵ100r̂k.
(A62)

Here, I=0.8227,ϵ100=r100/R0,r̂=r/r100,r̂k=rk/r100, and t̂=t/τA. Furthermore, δlinear(r) is the linear layer width (see Appendix 5), whereas ρ(r) is the plasma mass density profile.

As described in Ref. 24, the EPEC resonant plasma response model interpolates smoothly between the linear regime and the nonlinear Rutherford regime.27,53 The linear regime corresponds to Ŵkδ̂k, whereas the nonlinear region corresponds to Ŵkδ̂k.

5. Linear layer width

Let

τH(r)=R0B0gμ0ρns,
(A63)
τE(r)=r2D+(2/3)χe,
(A64)
τM(r)=R02q2χϕ,
(A65)
ρs(r)=miTeeB0g,
(A66)
τ(r)=ωiωe,
(A67)
S(r)=τRτH,
(A68)
PM(r)=τRτM,
(A69)
PE(r)=τRτE,
(A70)
D(r)=53S1/3ρsr,
(A71)
QE(r)=S1/3nωEτH,
(A72)
Qe,i(r)=S1/3nωe,iτH.
(A73)

Here, χe(r) and D(r) are the perpendicular electron energy and particle diffusivity profiles, respectively. According to the analysis of Ref. 52, the constant-ψ linear layer width is determined from the solution of

d2Ydp2[QE(QEQi)+i(QEQi)(PM+PE)p2+PMPEp4i(QEQe)+{PE+i(QEQi)D2}p2+(1+τ)PMD2p4]p2Y=0.
(A74)

If the small-p behavior of solution of the previous equation that is well-behaved as p is

Y(p)=Y0[1cp+O(p2)],
(A75)

then the linear layer width is

δlinearr=π|c|S1/3.
(A76)
6. Plasma angular velocity evolution

The quantity ϖ̂k that appears in Eqs. (A54) and (A55) evolves in time according to22 

ϖ̂k(t̂)=ϖ̂k0k=1,Kp=1,mkmkyp(r̂k)yp(r̂k)αk,p(t̂)k=1,Kp=1,zp(r̂k)zp(r̂k)βk,p(t̂).
(A77)

Here, ϖ̂k0=ϖk0τA,yp(r̂)=J1(j1,pr̂)/r̂, and zp(r̂)=J0(j0,pr̂). Moreover, ϖk0 is the so-called “natural frequency” (in the absence of the RMP) at the kth resonant surface; this quantity is defined as the helical phase velocity of a naturally unstable island chain, resonant at the surface, in the absence of an RMP (or any other island chains). Furthermore, Jm(z) is a standard Bessel function, and jm,p denotes the pth zero of this function. The time evolution equations for the αk,p and βk,p parameters specify how the plasma poloidal and toroidal angular velocity profiles are modified by the electromagnetic torques that develop within the plasma, in response to the applied RMP, and how these modifications affect the natural frequencies. The evolution equations take the form22 

(1+2Qk2)dαk,pdt̂+(j1,p2τ̂Mk+1τ̂θk+1τ̂cxk)αk,p=mk2[yp(r̂k)]2ρ̂kϵ1002[J2(j1,p)]2δT̂k,
(A78)
dβk,pdt̂+(j0,p2τ̂Mk+1τ̂cxk)βk,p=n2[zp(r̂k)]2ρ̂k[J1(j0,p)]2δT̂k,
(A79)

where

δT̂k=k=1,KÊkkΨ̂kΨ̂ksin(φkφkξkk)+ÊkkΨ̂kχ̂ksin(φkζk).
(A80)

Here, Qk=Q(rk),ρ̂k=ρ(rk)/ρ(0),τ̂Mk=r1002/[χϕ(rk)τA],τ̂θk=τθ(rk)/τA,τ̂cxk=τcx(rk)/τA. Moreover, χϕ(r) is the perpendicular toroidal momentum diffusivity profile, whereas

τθ(r)=τiiμ00i/(1+q2R02r2akk)
(A81)

is the poloidal flow damping timescale, and

τcx(r)=1nnσvicx
(A82)

is the charge exchange damping timescale. Furthermore,54 

Q2(r)=q2R022r2(1R21R2)/|r|2R2.
(A83)

Note that we are now using a more accurate expression for the enhancement of toroidal plasma inertia due to neoclassical poloidal flow damping, i.e., 1+2Qk2, rather than the large aspect-ratio approximation 1+2qk2. We have also incorporated charge exchange damping into the angular velocity evolution equations.

7. Natural frequencies

According to linear tearing mode theory, in the absence of the RMP, the natural frequency of the tearing mode resonant at the kth resonant surface is given by28,55–57

ϖlineark=n(ωE+ωe)rk.
(A84)

According to nonlinear tearing mode theory, in the absence of the RMP, the natural frequency of the tearing mode resonant at the kth resonant surface is given by22,58

ϖnonlineark=n(ωE+ωi+ωnci)rk,
(A85)

where

ωnci(r)=G00iiωE[L00iiL01ii(ηi1+ηi)]ωi[L00iIL01iI(ηI1+ηI)]ωI.
(A86)

Note the presence of the first term on the right-hand side of the previous equation, which was erroneously omitted in Ref. 22. This term relates to the influence of charge exchange with neutrals on the natural frequency.

Given that our resonant plasma response model interpolates smoothly between the linear regime (Ŵkδ̂k) and the nonlinear regime (Ŵkδ̂k), it makes sense to adopt a model for the natural frequency that does the same. Hence, the EPEC model for the natural frequency is

ϖk0=δ̂kϖlineark+Ŵkϖnonlinearkδ̂k+Ŵk.
(A87)

Note that this model differs from Eq. (29) of Ref. 22 because there is no regime, intermediate between the linear and nonlinear regimes, in which the natural frequency is determined solely by the local E×B velocity. The reason for making this change is the absence of a solid theoretical basis for supposing that such a regime exists (unlike the linear and nonlinear regimes).

8. Island widths and temperature and density flattening widths

The full width (in ΨN) of the magnetic separatrix of the island chain induced at the kth resonant surface is

Wk=4[Ψ̂kσk|Ψ̂p(0)|]1/2,
(A88)

where σk=σ(rk),σ(r)=dlnq/dΨN, and Ψ̂p(r)=Ψp/(R02B0). However, the halfway point between the inner and outer extremities of the separatrix is located at the flux-surface where ΨN=ΨN(rk)AkWk2/8 and Ak=[(1/3)(d2q/dΨN2)/(dq/dΨN)]rk. On average, the presence of a magnetic island chain of full width Wk at the kth resonant surface causes the electron temperature, ion temperature, and electron number density profiles to be flattened in annular regions, centered on the halfway points between the inner and outer radii of the island separatrices, of widths (in ΨN),22,59–61

δTek=2πWktanh(WkWTek),
(A89)
δTik=2πWktanh(WkWTik),
(A90)
δnek=2πWktanh(WkWnek),
(A91)

respectively, where

WTek=(χevTeR0f2nσ|Ψ̂p(0)|2)rk1/3,
(A92)
WTik=(χivTiR0f2nσ|Ψ̂p(0)|2)rk1/3,
(A93)
Wnek=(DvTiR0f2nσ|Ψ̂p(0)|2)rk1/3.
(A94)

Here, χi(r) is the perpendicular ion energy diffusivity profile. Note that the previous two expression take into account the fact that the parallel transport of both energy and particles in the pedestal region of a high temperature tokamak plasma is convective, rather than diffusive, in nature.59,60

The total change in the plasma pressure interior to the innermost resonant surface due to the temperature and density flattening at the various resonant surfaces (assuming that the pressure at the LCFS is fixed) is61 

ΔP=k=1,K{[δnekdnedΨNTeδnekWk2Ak8(dnedΨNdTedΨN+d2nedΨn2Te)+δnek3d3nedΨN3Te24]+[δTeknedTedΨNδTekWk2Ak8(dnedΨNdTedΨN+ned2TedΨn2)+δTek3ne24d3TedΨN3]+[δnekdnidΨNTiδnekWk2Ak8(dnidΨNdTidΨN+d2nidΨn2Ti)+δnek3d3nidΨN3Ti24]+[δTiknidTidΨNδTikWk2Ak8(dnidΨNdTidΨN+nid2TidΨn2)+δTik3ni24d3TidΨN3]}rk,
(A95)

where we have neglected the impurity ion pressure.

1.
F.
Wagner
,
G.
Becker
,
K.
Behringer
,
D.
Campbell
,
A.
Eberhagen
,
W.
Engelhardt
,
G.
Fussmann
,
O.
Gehre
,
J.
Gernhardt
,
G. V.
Gierke
 et al,
Phys. Rev. Lett.
49
,
1408
(
1982
).
2.
H.
Zohm
,
Plasma Phys. Controlled Fusion
38
,
105
(
1996
).
3.
A.
Loarte
,
G.
Saibene
,
R.
Sartori
,
M.
Bécoulet
,
L.
Horton
,
T.
Eich
,
A.
Herrmann
,
M.
Laux
,
G.
Matthews
,
S.
Jachmich
 et al,
J. Nucl. Mater.
313–316
,
962
(
2003
).
4.
T. E.
Evans
,
R. A.
Moyer
,
J. G.
Watkins
,
P. R.
Thomas
,
T. H.
Osborne
,
J. A.
Boedo
,
M. E.
Fenstermacher
,
K. H.
Finken
,
R. J.
Groebner
,
M.
Groth
 et al,
Phys. Rev. Lett.
92
,
235003
(
2004
).
5.
Y.
Liang
,
H. R.
Koslowski
,
P. R.
Thomas
,
E.
Nardon
,
B.
Alper
,
P.
Andrew
,
Y.
Andrew
,
G.
Arnoux
,
Y.
Baranov
,
M.
Bécoulet
 et al,
Phys. Rev. Lett.
98
,
265004
(
2007
).
6.
W.
Suttrop
,
T.
Eich
,
J. C.
Fuchs
,
S.
Günter
,
A.
Janzer
,
A.
Herrmann
,
A.
Kallenbach
,
P. T.
Lang
,
T.
Lunt
,
M.
Maraschek
 et al,
Phys. Rev. Lett.
106
,
225004
(
2011
).
7.
Y. M.
Jeon
,
J.-K.
Park
,
S. W.
Yoon
,
W. H.
Ko
,
S. G.
Lee
,
K. D.
Lee
,
G. S.
Yun
,
Y. U.
Nam
,
W. C.
Kim
,
J.-G.
Kwak
,
K. S.
Lee
,
H. K.
Kim
,
H. L.
Yang
 et al,
Phys. Rev. Lett.
109
,
035004
(
2012
).
8.
A.
Kirk
,
I. T.
Chapman
,
Y.
Liu
,
P.
Cahyna
,
P.
Denner
,
G.
Fishpool
,
C. J.
Ham
,
J. R.
Harrison
,
Y.
Liang
,
E.
Nardon
,
S.
Saarelma
,
R.
Scannell
,
A. J.
Thornton
, and
MAST Team
,
Nucl. Fusion
53
,
043007
(
2013
).
9.
Y.
Sun
,
Y.
Liang
,
Y. Q.
Liu
,
S.
Gu
,
X.
Yang
,
W.
Guo
,
T.
Shi
,
M.
Jia
,
L.
Wang
,
B.
Lyu
 et al,
Phys. Rev. Lett.
117
,
115001
(
2016
).
10.
T. E.
Evans
,
M. E.
Fenstermacher
,
R. A.
Moyer
,
T. H.
Osbourne
,
J. G.
Watkins
,
P.
Gohil
,
I.
Joseph
,
M. J.
Schaffer
,
L. R.
Baylor
 et al,
Nucl. Fusion
48
,
024002
(
2008
).
11.
M.
Kim
,
J.
Lee
,
W. H.
Ho
,
S.-H.
Hahn
,
Y.
In
,
Y. M.
Jeon
,
W.
Suttrop
,
S. K.
Kim
,
G. Y.
Park
,
J.-W.
Juhn
, and
J. H.
Lee
,
Phys. Plasmas
27
,
112501
(
2020
).
12.
Q. M.
Hu
,
R.
Nazikian
,
B. A.
Grierson
,
N. C.
Logan
,
J.-K.
Park
,
C.
Paz-Solden
, and
Q.
Yu
,
Phys. Plasmas
26
,
120702
(
2019
).
13.
Q. M.
Hu
,
R.
Nazikian
,
B. A.
Grierson
,
N. C.
Logan
,
C.
Paz-Solden
, and
Q.
Yu
,
Nucl. Fusion
60
,
076001
(
2020
).
14.
Q. M.
Hu
,
R.
Nazikian
,
B. A.
Grierson
,
N. C.
Logan
,
D. M.
Orlov
,
C.
Paz-Solden
, and
Q.
Yu
,
Phys. Rev. Lett.
125
,
045001
(
2020
).
15.
Q. M.
Hu
,
R.
Nazikian
,
N. C.
Logan
,
J.-K.
Park
,
C.
Paz-Solden
,
S. M.
Yang
,
B. A.
Grierson
,
Y.
In
,
Y. M.
Jeon
,
M.
Kim
,
S. K.
Kim
,
D. M.
Orlov
,
G. Y.
Park
, and
Q.
Yu
,
Phys. Plasmas
28
,
052505
(
2021
).
16.
Q.
Yu
,
S.
Günter
, and
B. D.
Scott
,
Phys. Plasmas
10
,
797
(
2003
).
17.
Q.
Yu
,
Nucl. Fusion
50
,
025014
(
2010
).
18.
Q.
Yu
and
S.
Günter
,
Nucl. Fusion
51
,
073030
(
2011
).
19.
J.-K.
Park
and
N. C.
Logan
,
Phys. Plasmas
24
,
032505
(
2017
).
20.
P. B.
Snyder
,
T. H.
Osbourne
,
K. H.
Burrell
,
R. J.
Groebner
,
A. W.
Leonard
,
R.
Nazikian
,
D. M.
Orlov
,
O.
Schmitz
,
M. R.
Wade
, and
H. R.
Wilson
,
Phys. Plasmas
19
,
056115
(
2012
).
21.
R.
Fitzpatrick
,
Nucl. Fusion
33
,
1049
(
1993
).
22.
R.
Fitzpatrick
and
A. O.
Nelson
,
Phys. Plasmas
27
,
072501
(
2020
).
23.
R.
Fitzpatrick
,
Phys. Plasmas
27
,
102511
(
2020
).
24.
R.
Fitzpatrick
,
Phys. Plasmas
28
,
022503
(
2021
).
25.
H. P.
Furth
,
J.
Killeen
, and
M. N.
Rosenbluth
,
Phys. Fluids
6
,
459
(
1963
).
26.
B.
Coppi
,
J. M.
Greene
, and
J. L.
Johnson
,
Nucl. Fusion
6
,
101
(
1966
).
27.
P. H.
Rutherford
,
Phys. Fluids
16
,
1903
(
1973
).
28.
G.
Ara
,
B.
Basu
,
B.
Coppi
,
G.
Laval
,
M. N.
Rosenbluth
, and
B. V.
Waddell
,
Ann. Phys.
112
,
443
(
1978
).
29.
A.
Pletzer
and
R. L.
Dewar
,
J. Plasma Phys.
45
,
427
(
1991
).
30.
R.
Fitzpatrick
,
R. J.
Hastie
,
T. J.
Martin
, and
C. M.
Roach
,
Nucl. Fusion
33
,
1533
(
1993
).
31.
S.
Tokuda
,
Nucl. Fusion
41
,
1037
(
2001
).
32.
D. P.
Brennan
,
R. J.
La Haye
,
A. D.
Turnbull
,
M. S.
Chu
,
T. H.
Jensen
,
L. L.
Lao
,
T. C.
Luce
,
P. A.
Politzer
, and
E. J.
Strait
,
Phys. Plasmas
10
,
1643
(
2003
).
33.
A. H.
Glasser
,
Z. R.
Wang
, and
J.-K.
Park
,
Phys. Plasmas
23
,
112506
(
2016
).
34.
R.
Fitzpatrick
,
Phys. Plasmas
24
,
072506
(
2017
).
35.
H. K.
Kim
,
H. L.
Yang
,
G. H.
Kim
,
J. Y.
Kim
,
H.
Jhang
,
J. S.
Bak
, and
G. S.
Lee
,
Fusion Eng. Des.
84
,
1029
(
2009
).
36.
J.
Lee
,
Y. M.
Jeon
,
Y.
In
,
G. Y.
Park
,
G. S.
Yun
,
W.
Lee
,
M.
Kim
,
J. H.
Lee
,
W. H.
Ko
,
H. K.
Park
 et al,
Nucl. Fusion
59
,
066033
(
2019
).
37.
S. K.
Kim
,
S.
Pamela
,
O.
Kwon
,
M.
Bécoulet
,
G. T. A.
Huysmans
,
Y.
In
,
M.
Hoelzl
,
J. H.
Lee
,
M.
Kim
,
G. Y.
Park
 et al,
Nucl. Fusion
60
,
026009
(
2020
).
38.
P.
Monier-Garbet
,
K. H.
Burrell
,
F. L.
Hinton
,
J.
Kim
,
X.
Garbet
, and
R. J.
Groebner
,
Nucl. Fusion
37
,
403
(
1997
).
39.
C. F.
Barnett
, “
Atomic data for fusion
,”
Report No. ORNL-6086/V1
(
Oak Ridge National Laboratory
,
TN
,
1990
).
40.
S. P.
Hirshman
and
D. J.
Sigmar
,
Nucl. Fusion
21
,
1079
(
1981
).
41.
H.
Lütjens
,
A.
Bondeson
, and
O.
Sauter
,
Comp. Phys. Commun.
97
,
219
(
1996
).
42.
J. W.
Connor
,
R. J.
Hastie
,
H. R.
Wilson
, and
R. L.
Miller
,
Phys. Plasmas
5
,
2687
(
1998
).
43.
P. B.
Snyder
,
H. R.
Wilson
,
J. R.
Ferron
,
D.
Mossessian
,
M.
Murakami
, and
X. Q.
Xu
,
Phys. Plasmas
9
,
2037
(
2002
).
44.
Y.
In
,
A.
Loarte
,
H. H.
Lee
,
K.
Kim
,
Y. M.
Jeon
,
J.-K.
Park
,
J.-W.
Ahn
,
G. Y.
Park
,
M.
Kim
,
H.
Park
 et al,
Nucl. Fusion
59
,
126045
(
2019
).
45.
Y.
In
,
Y. M.
Jeon
,
J.-K.
Park
,
A.
Loarte
,
J.-W.
Ahn
,
J. H.
Lee
,
H. H.
Lee
,
G. Y.
Park
 et al,
Nucl. Fusion
59
,
056009
(
2019
).
46.
Y.
In
,
J.-K.
Park
,
Y. M.
Jeon
,
J.
Kim
,
G. Y.
Park
,
J.-W.
Ahn
,
A.
Loarte
,
W. H.
Ko
 et al,
Nucl. Fusion
57
,
116054
(
2017
).
47.
O.
Meneghini
,
S. P.
Smith
,
L. L.
Lao
,
O.
Izacard
,
Q.
Ren
,
J. M.
Park
,
J.
Candy
 et al,
Nucl. Fusion
55
,
083008
(
2015
).
48.
J.-K.
Park
,
M. J.
Schaffer
,
J. E.
Menard
, and
A. H.
Boozer
,
Phys. Rev. Lett.
99
,
195003
(
2007
).
49.
A. H.
Boozer
and
C.
Nührenberg
,
Phys. Plasmas
13
,
102501
(
2006
).
50.
Y. B.
Kim
,
P. H.
Diamond
, and
R. J.
Groebner
,
Phys. Fluids B
3
,
2050
(
1991
).
51.
Y. B.
Kim
,
P. H.
Diamond
, and
R. J.
Groebner
,
Phys. Fluids B
4
,
2996
(
1992
).
52.
A.
Cole
and
R.
Fitzpatrick
,
Phys. Plasmas
13
,
032503
(
2006
).
53.
P. H.
Rutherford
, “
Basic physical processes of toroidal fusion plasmas
,” in
Proceedings of Course and Workshop
, Varenna (
Commission of the European Communities
,
Brussels
,
1985
), Vol.
2
, p.
531
.
54.
S. P.
Hirshman
,
Nucl. Fusion
18
(
7
),
917
(
1978
).
55.
M.
Bécoulet
,
F.
Orain
,
P.
Maget
,
N.
Mellet
,
X.
Garbet
,
E.
Nardon
,
G. T. A.
Huysmans
,
T.
Caspar
,
A.
Loarte
,
P.
Cayna
 et al,
Nucl. Fusion
52
,
054003
(
2012
).
56.
N. M.
Ferraro
,
Phys. Plasmas
19
,
056105
(
2012
).
57.
F.
Orain
,
M.
Bécoulet
,
G.
Dif-Pradalier
,
G. T. A.
Huysmans
,
S.
Pamela
,
E.
Nardon
,
C.
Passeron
,
G.
Latu
,
V.
Grandgirard
,
A.
Fil
 et al,
Phys. Plasmas
20
,
102510
(
2013
).
58.
R.
Fitzpatrick
,
Phys. Plasmas
25
,
082513
(
2018
).
59.
R.
Fitzpatrick
,
Phys. Plasmas
2
,
825
(
1995
).
60.
N. N.
Gorelenkov
,
R. V.
Budny
,
Z.
Chang
,
M. V.
Gorelenkova
, and
L. E.
Zakharov
,
Phys. Plasmas
3
,
3379
(
1996
).
61.
Z.
Chang
and
J. D.
Callen
,
Nucl. Fusion
30
,
219
(
1990
).