The EPEC code is employed to model the *q*_{95} 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 *q*_{95} (i.e., $\Delta q95\u22430.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)].

## I. INTRODUCTION

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 EAST^{9} tokamaks.

RMP-induced ELM suppression is only observed to occur when *q*_{95} [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., $\Delta q95\u22430.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 *q*_{95}. Recent simulations^{12–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 *q*_{95} 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 region^{20} (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\varphi /B\theta )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 damping^{21}) 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 code^{22–24} employs an asymptotic matching^{21,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 *q*_{95} 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.

## II. KSTAR DISCHARGE #18594

### A. Introduction

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.

### B. Plasma equilibrium

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=\u22121.79\u2009T,\u2009R0=1.80\u2009m,\u2009r100=0.598\u2009m$, $q95=4.04$ (these quantities are defined in Appendix 2 a and 2 b), and ELM suppression has not yet occurred.

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 $1\u2009m2\u2009s\u22121$, whereas the perpendicular particle diffusivity has been given the plausible value $0.20\u2009m2\u2009s\u22121$. It can be seen that the top of the pedestal corresponds to $\Psi N=0.945$. (The normalized poloidal flux $\Psi N$ is defined in Appendix 2 b.)

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

(The flux-surface label *r* is defined in Appendix 2 a. *r*_{100} is the flux-surface label of the LCFS. The flux-surface average operator $\u27e8\cdots \u27e9$ is defined in Appendix 3 e.) We shall consider two cases: a “high-recycling” case in which $\u27e8nn\u27e9(r100)=1.3\xd71017\u2009m\u22123$ (which is similar to the inferred DIII-D value), and a “low-recycling” case in which $\u27e8nn\u27e9(r100)=2.5\xd71016\u2009m\u22123$. The parameter *l _{n}* is given the value $ln=0.013\u2009m$.

^{38}The flux-surface neutral poloidal asymmetry parameter (see Appendix 3 e),

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 value^{39}

### C. E × B frequency profile

Unfortunately, there are no usable poloidal rotation data in KSTAR discharge #18594. Hence, it is necessary to deduce the $E\xd7B$ 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\xd7B$ frequency profile is related to the impurity ion toroidal angular velocity profile as follows (see Appendix 3 f),

Here, $\omega \varphi \u200aI,\u2009\omega \theta \u200aI$, and $\omega \u2217\u200aI$ 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\xd7B$ frequency. In the low-recycling case, there is an

*E*well in the pedestal (where the

_{r}*ω*curve dips below zero). However, in the high-recycling case, this well is eliminated via the action of charge exchange with neutrals.

_{E}### D. Natural frequencies

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.

## III. RESCALING OF PLASMA EQUILIBRIUM

In order to search for *q*_{95} windows in which RMP-induced ELM suppression is possible, we require a collection of otherwise similar plasma equilibria with a range of different *q*_{95} 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

where

Here, $R0\u2009B0\u2009g(r)\u2009\u2207\varphi $ is the toroidal magnetic field-strength, where $\varphi $ is a geometric toroidal angle (see Appendix 2 b). Moreover, *r*_{95} is the value of the flux-surface label *r* at the surface at which $\Psi N=0.95$ (see Appendix 2 b). The second rescaling is such that

where

Here, the species label *a* stands for either *e* (electrons), *i* (majority ions), or *I* (impurity ions) (see Appendix 3 a). Moreover, $\Psi p(r)$ is the poloidal magnetic flux (divided by $2\pi $) (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 *q*_{95} is changed from its original value of 4.04 to the new value $q95\u200atarget$, 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 model^{43} into the rescaling method, in order to avoid this possibility. This will be attempted in future publications.

## IV. OPTIMUM RMP COIL PHASING

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

Here, $\Delta \Psi 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, *E _{kk}* is the (real) normalized tearing stability index at the resonant surface,

^{25}and

*χ*is a (complex) parameter determined by the GPEC code (see Appendix 2 d) that measures the ideal response of the

_{k}*k*th resonant surface to the currents flowing in the RMP coils. In fact, the previous equation can also be written

Here, $Irmp$ is the (real) peak current flowing in the upper, middle, and lower coil sets (in units of kiloamps per turn). Moreover, $\chi kU,\u2009\chi kM$, and $\chi 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, $\Delta \varphi $ 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

*n*th toroidal harmonic of the current pattern flowing in the middle coil set varies with toroidal angle as $cos\u2009(n\u2009\varphi )$, then that in the upper coil set varies as $cos\u2009[n\u2009(\varphi \u2212\Delta \varphi )]$, and that in the lower coil set as $cos\u2009[n\u2009(\varphi +\Delta \varphi )]$.} 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., $\Psi N=0.945$). Suppose that the closest resonant surface to the top of the pedestal is the *k*th 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 $|\Delta \Psi k|$. (In fact, we maximize the modulus of a linear superposition of the $\Delta \Psi 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

where $\gamma k\u200aUM,\u2009\gamma k\u200aUL$, and $\gamma k\u200aML$ are real, then the optimization procedure boils down to maximizing the function

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 *q*_{95}. (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\xb0$, which is in accordance with KSTAR experimental observations.^{37} However, the actual optimal phase-shift is slightly less than $90\xb0$ for large *q*_{95} (i.e., $q95>3.7$), and slightly larger for small *q*_{95}. The optimal phase-shift for *n* = 1 RMPs exhibits a much greater variation with *q*_{95}. The phase-shift used in experiments is $270\xb0$.^{44} It can be seen that the optimal phase-shift determined by EPEC approaches this value at small *q*_{95}, but is significantly different at large *q*_{95}.

## V. DETERMINATION OF *q*_{95} WINDOWS FOR RMP-INDUCED ELM SUPPRESSION

### A. *n* = 2

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=4\u2009kA/turn$ is applied to each row of RMP coils with the optimal phasing. Using the rescaled plasma equilibria described in Sec. III, *q*_{95} 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<\Psi N<1.00$ are neglected in the calculation (because GPEC does not give reliable results for such surfaces).

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 *q*_{95}, 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 *q*_{95} 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 $5\u2009krad/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.

Figures 8 and 9 show the full widths (in $\Psi 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 *q*_{95}, 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 ($\Delta \Psi N\u22430.1$) magnetic island chain at the top of the pedestal. Again, the differences between the low- and high-recycling cases is comparatively minor.

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\u2032$ 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., $\Psi N=0.945$), calculated as functions of *q*_{95} 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 *q*_{95} lies in certain narrow (i.e., $\Delta q95\u22430.1$) windows that are evenly spaced in *q*_{95} 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 *q*_{95}.

Now, the experimental $m=8/n=2$ RMP-induced ELM suppression window in KSTAR occurs at $q95\u22433.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, *q*_{95} 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 *q*_{95} 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 *q*_{95} 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 *q*_{95} windows and experimental observations. (It should be noted that the *q*_{95} windows in Fig. 8 of Ref. 15 have been shifted down by 0.2 in order to take into account the difference between the *q*_{95} 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\xd7B$ velocity profile in the pedestal, which, in turn, leads to a shift in the locations of the *q*_{95} 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 *q*_{95}) 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

*E*, 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\u2032$. 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.

_{r}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 *q*_{95} 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 *q*_{95}. 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 *q*_{95} 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 *q*_{95}), widths (in *q*_{95}), 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 *q*_{95} 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.

### B. *n* = 1

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, *q*_{95} is ramped downward from a value of about 5.5 to a value of about 2.5 over a time interval of 3 s.

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 *q*_{95} lies in certain narrow (i.e., $\Delta q95\u22430.1$) windows that are evenly spaced in *q*_{95} 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 *q*_{95}.

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 *q*_{95} 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 *q*_{95}.

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 *q*_{95} 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 *q*_{95}) 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 *q*_{95}) 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}

## VI. SUMMARY AND CONCLUSIONS

We have employed the EPEC code (see the Appendix) to model the *q*_{95} 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 *q*_{95} vs $Irmp$ space. Here, *q*_{95} 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 *q*_{95} (i.e., $\Delta q95\u22430.1$) and that (assuming a relative pedestal pressure reduction of 15% is required to trigger ELM suppression^{15}) 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 *q*_{95} 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 *q*_{95}), height (in $Irmp$), and position (in *q*_{95})] 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\xd7B$ velocity profile in the pedestal, which, in turn, leads to a shift in the locations of the *q*_{95} windows. Unfortunately, data on neutral particles within the LCFS of H-mode tokamak plasmas are presently very hard to come by.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

### APPENDIX: DESCRIPTION OF EPEC MODEL

##### 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*, $\varphi $, and *Z* be right-handed cylindrical coordinates whose symmetry axis corresponds to the toroidal symmetry axis of the plasma. Let *r*, *θ*, and $\varphi $ be right-handed flux coordinates whose Jacobian is $J\u2261(\u2207r\xd7\u2207\theta \xb7\u2207\varphi )\u22121=r\u2009R\u200a2/R0$. Here, *R*_{0} 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\pi $ 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<\theta <\pi $ correspond to the region above the midplane.

###### b. Equilibrium magnetic field

The equilibrium magnetic field is written $B=R0\u2009B0[f(r)\u2009\u2207\varphi \xd7\u2207r+g(r)\u2009\u2207\varphi ]$, where *B*_{0} is a convenient scale toroidal magnetic field-strength, and $q(r)=r\u2009g/(R0\u2009f)$ is the safety-factor profile. The equilibrium poloidal magnetic flux (divided by $2\pi $), $\Psi p(r)$, satisfies $d\Psi p/dr=R0\u2009B0\u2009f(r)$, where, by convention, $\Psi p(r100)=0$. The normalized poloidal magnetic flux, $\Psi N(r)$, is defined such that $\Psi N(r)=1\u2212\Psi p(r)/\Psi p(0)$. Hence, $\Psi N(0)=0$ and $\Psi N(r100)=1$. Finally, if $\Psi N(r95)=0.95$, then $q95\u2261q(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 form^{22}

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 *k*th resonant surface, *r* = *r _{k}*, at which $n\u2009q(rk)=mk$, where

*m*is a positive integer. Let $\Psi k=\psi k(rk)/mk$, and $\Delta \Psi k=[\Xi k]rk\u2212rk+$.

_{k}^{22}Here, $\Psi k$ is the (complex) reconnected helical magnetic flux (divided by $2\pi \u2009R0$) at the

*k*th resonant surface, whereas $\Delta \Psi k$ (which has the same dimensions as $\Psi 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 $\Psi k$ and the $\Delta \Psi k$ values are related according to the inhomogeneous toroidal tearing mode dispersion relation, which takes the form^{21,22,29–34}

Here, $Ekk\u2032$ (for *k*, $k\u2032=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\u2032$ matrix using a high-*q* approximation. In fact, if $Fkk\u2032$ is the inverse of the $Ekk\u2032$ matrix, then^{22}

and

with

Here, the double integral in Eq. (A5) is taken around the *k*th resonant surface (cylindrical coordinates *R _{k}*, 0,

*Z*; flux coordinates

_{k}*r*,

_{k}*θ*, 0, with

_{k}*r*constant; resonant poloidal mode number

_{k}*m*) and the $k\u2032$ th resonant surface (cylindrical coordinates $Rk\u2032$, 0, $Zk\u2032$; flux coordinates $rk\u2032,\u2009\theta k\u2032$, 0, with $rk\u2032$ constant; resonant poloidal mode number $mk\u2032$). Finally, the $\Gamma (z)$ and $P\mu \u200a\nu (z)$ are the gamma functions and associated Legendre functions, respectively.

_{k}The (complex) *χ _{k}* parameters are determined from the GPEC code.

^{19,48}To be more exact, the GPEC code calculates the (complex) dimensionless $\Delta mk\u200an$ 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 $\Delta mk\u200an$ parameters are related to the

*χ*parameters according to

_{k}^{22}

where $qk=mk/n$, and $akk(r)=\u222e|\u2207r|\u22122\u2009d\theta /(2\pi )$. 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=\u2212e$, *e _{i}* =

*e*, and $eI=ZI\u2009e$, respectively, where

*e*is the magnitude of the electron charge. Quasi-neutrality demands that $ne=ni+ZI\u2009nI$, where $na(r)$ is the species-

*a*number density. Let $\alpha I(r)=ZI\u2009(Zeff\u22121)/(ZI\u2212Zeff)$, where $Zeff(r)=(ni+ZI\u200a2\u2009nI)/ne$ is the effective ion charge number. It follows that $ni/ne=(ZI\u2212Zeff)/(ZI\u22121)$ and $nI/ne=(Zeff\u22121)/[ZI\u2009(ZI\u22121)]$.

###### b. Collisionality parameters

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

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

Let

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 $\omega t\u200aa(r)=Kt\u2009\gamma \u2009vT\u200aa$, where

and $vT\u200aa=2\u2009Ta/ma$.^{40} Here, *m _{a}* is the species-

*a*mass, and $Ta(r)$ the species-

*a*temperature (in energy units). The fraction of circulating particles is

^{40}

Finally, the dimensionless species-*a* collisionality parameter is written^{40} $\nu \u2217\u200aa(r)=K\u2217\u2009gt/(\omega t\u200aa\u2009\tau aa)$, where $gt(r)=fc/(1\u2212fc)$,

Here, the Coulomb logarithm, $ln\u2009\Lambda $, is assumed to take the same large constant value (i.e., $ln\u2009\Lambda \u224317$), independent of the species.

###### c. Collisional friction matrices

Let $xab=vT\u200ab/vT\u200aa$. In the following, all quantities that are of order $(me/mi)\u200a1/2,\u2009(me/mI)\u200a1/2$, or smaller, are neglected with respect to unity. The 2 × 2 dimensionless ion collisional friction matrices, $[F\u200aii](r),\u2009[F\u200aiI](r),\u2009[F\u200aIi](r)$, and $[F\u200aII](r)$, are defined to have the following elements:^{40,50,51}

The 2 × 2 dimensionless electron collisional friction matrix, $[F\u200aee](r)$, is defined to have the following elements:^{40,50,51} $F00\u200aee=Zeff$, $F01\u200aee=(3/2)\u2009Zeff$, $F10\u200aee=F01\u200aee$, $F11\u200aee=2+(13/4)\u2009Zeff$.

###### d. Neoclassical viscosity matrices

The 2 × 2 dimensionless species-*a* neoclassical viscosity matrix, $[\mu \u200aa](r)$, is defined to have the following elements: $\mu 00\u200aa=K00\u200aa,\u2009\mu 01\u200aa=(5/2)\u2009K00\u200aa\u2212K01\u200aa$, $\mu 10\u200aa=\mu 01\u200aa$, $\mu 11\u200aa=K11\u200aa\u22125\u2009K01\u200aa+(25/4)\u2009K00\u200aa$.^{40,50,51} Here,

and

Furthermore,

and, finally,

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 $[\mu \u0303\u200aI]=\alpha I\u200a2\u2009(Ti/TI)\u2009xIi\u2009[\mu \u200aI]$. The requirement of equilibrium force and heat balance parallel to the magnetic field^{38,40,50,51} leads us to define four 2 × 2 dimensionless ion matrices, $[L\u200aii](r),\u2009[L\u200aiI](r),[L\u200aIi](r)$, and $[L\u200aII](r)$, where

the additional four 2 × 2 dimensionless ion matrices, $[G\u200aii](r),\u2009[G\u200aiI](r),\u2009[G\u200aIi](r)$, and $[G\u200aII](r)$, where

and the 2 × 2 dimensionless electron matrix, $[Q\u200aee](r)$, where $[Q\u200aee]=[F\u200aee+\mu \u200ae]\u200a\u22121$. Here,^{38}

where

Moreover, $\u27e8\sigma \u2009v\u27e9i\u200acx$ is the flux-surface averaged rate constant for charge exchange reactions between neutrals and majority ions, $nn(r,\Theta )$ the neutral particle number density, and $En/Ti$ the ratio of the incoming neutral energy to the majority ion energy. The parameter *y _{n}* 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

where $V\u200aI$ is the impurity ion fluid velocity, and the right-hand sides are evaluated on the outboard midplane. According to neoclassical theory^{40}

where

and

Here, $pa(r)=na\u2009Ta$, and $\Phi (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 $\Psi k=R0\u2009B0\u2009\Psi \u0302k\u2009e\u2212i\u200a\phi k,\u2009\chi k=R0\u2009B0\u2009\chi \u0302k\u2009e\u2212i\u200a\zeta k$, and $Ekk\u2032=E\u0302kk\u2032\u2009e\u2212i\u200a\xi kk\u2032$, where $\Psi \u0302k>0,\u2009\phi k,\u2009\chi \u0302k>0$, *ζ _{k}*, $E\u0302kk\u2032>0$, and $\xi kk\u2032$ are all real quantities. Furthermore, let $Xk=\Psi \u0302k\u2009\u2009cos\u2009\phi k$ and $Yk=\Psi \u0302k\u2009\u2009sin\u2009\phi k$. The resonant plasma response model at the

*k*th resonant surface takes the form

^{24}

where

Here, $I=0.8227,\u2009\u03f5100=r100/R0,\u2009r\u0302=r/r100,\u2009r\u0302k=rk/r100$, and $t\u0302=t/\tau A$. Furthermore, $\delta linear(r)$ is the linear layer width (see Appendix 5), whereas $\rho (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 $W\u0302k\u226a\delta \u0302k$, whereas the nonlinear region corresponds to $W\u0302k\u226b\delta \u0302k$.

##### 5. Linear layer width

Let

Here, $\chi e(r)$ and $D\u22a5(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

If the small-*p* behavior of solution of the previous equation that is well-behaved as $p\u2192\u221e$ is

then the linear layer width is

##### 6. Plasma angular velocity evolution

The quantity $\varpi \u0302k$ that appears in Eqs. (A54) and (A55) evolves in time according to^{22}

Here, $\varpi \u0302k\u200a0=\varpi k\u200a0\u2009\tau A,\u2009yp(r\u0302)=J1(j1,p\u2009r\u0302)/r\u0302$, and $zp(r\u0302)=J0(j0,p\u2009r\u0302)$. Moreover, $\varpi k\u200a0$ is the so-called “natural frequency” (in the absence of the RMP) at the *k*th 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 *p*th zero of this function. The time evolution equations for the $\alpha k,p$ and $\beta 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 form^{22}

where

Here, $Qk=Q(rk),\u2009\rho \u0302k=\rho (rk)/\rho (0),\u2009\tau \u0302M\u200ak=r100\u200a2/[\chi \varphi (rk)\u2009\tau A],\u2009\tau \u0302\theta \u200ak=\tau \theta (rk)/\tau A,\u2009\tau \u0302cx\u200ak=\tau cx(rk)/\tau A$. Moreover, $\chi \varphi (r)$ is the perpendicular toroidal momentum diffusivity profile, whereas

is the poloidal flow damping timescale, and

is the charge exchange damping timescale. Furthermore,^{54}

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+2\u2009Qk\u200a2$, rather than the large aspect-ratio approximation $1+2\u2009qk\u200a2$. 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 *k*th resonant surface is given by^{28,55–57}

According to nonlinear tearing mode theory, in the absence of the RMP, the natural frequency of the tearing mode resonant at the *k*th resonant surface is given by^{22,58}

where

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 ($W\u0302k\u226a\delta \u0302k$) and the nonlinear regime ($W\u0302k\u226b\delta \u0302k$), it makes sense to adopt a model for the natural frequency that does the same. Hence, the EPEC model for the natural frequency is

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\xd7B$ 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 $\Psi N$) of the magnetic separatrix of the island chain induced at the *k*th resonant surface is

where $\sigma k=\sigma (rk),\u2009\sigma (r)=d\u2009ln\u2009q/d\Psi N$, and $\Psi \u0302p(r)=\Psi p/(R0\u200a2\u2009B0)$. However, the halfway point between the inner and outer extremities of the separatrix is located at the flux-surface where $\Psi N=\Psi N(rk)\u2212Ak\u2009Wk\u200a2/8$ and $Ak=[(1/3)\u2009(d\u200a2q/d\Psi N\u200a2)/(dq/d\Psi N)]rk$. On average, the presence of a magnetic island chain of full width *W _{k}* at the

*k*th 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 $\Psi N$),

^{22,59–61}

respectively, where

Here, $\chi 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) is^{61}

where we have neglected the impurity ion pressure.