Transition in particle transport under resonant magnetic perturbations in a tokamak

Nonlinear 3D MHD simulations and validations reveal that the hybrid particle-MHD transport is a key process for driving the pump-out in the presence of Resonant Magnetic Perturbations (RMPs) in the KSTAR tokamak. Particle transport and the resulting density pump-out by RMPs are shown to be composed of not only the classical flow convection near magnetic islands due to polarization but also the neoclassical ion diffusion across perturbed magnetic surfaces. The latter is known as the Neoclassical Toroidal Viscosity (NTV) and is integrated into nonlinear MHD simulations here for the first time, revealing that the two-stage pump-outs observed in KSTAR experiments are reproduced only with such integrated nonlinear MHD and transport evolution. Near-resonant responses, which have received less attention than the resonant response, play distinct roles in the pump-out along with the island formation. In addition, this modeling is used to investigate the pump-outs in double-null-like plasmas and numerically capture the effect of the double-null shape on the pump-outs, which may explain the difficulty of Edge Localized Mode (ELM) suppression access in double-like plasmas. This reveals new aspects of the impact toroidal geometry and mode coupling have on 3D physics and reveals the importance of near-resonant components in suppressing ELMs.


Introduction
When sufficiently heated, magnetically confined tokamak plasmas spontaneously access to a high confinement mode (Hmode) [1]-a promising plasma operation scenario for future fusion power plants. The H-mode is characterized by a narrow edge transport barrier concomitant with the formation of an edge pedestal with a steep pressure gradient. This pedestal greatly improves global plasma confinement, but it can also lead to the growth of dangerous edge instabilities called Edge Localized Modes (ELMs) [2]. ELMs trigger rapid relaxation of the edge density and temperature, resulting in intense transient heat fluxes on the reactor walls and undesired material erosion [3]. Therefore, in order to retain a tokamak design as a viable fusion reactor, it is crucial to find ways to achieve high plasma confinement while simultaneously suppressing ELM events.
One of the most effective methods to control ELMs is to apply Resonant Magnetic Perturbations (RMPs) using 3D coils [4]. RMPs purposefully cause additional transport in the pedestal, degrading its height to a point where ELMs are fully suppressed [5][6][7][8][9][10][11]. Notably, applying RMP leads to a substantial reduction of density pedestal called pumpout, considered a necessary evil for accessing RMP-induced ELM-free plasmas. Extensive studies have been conducted to understand and predict the RMP pump-out phenomena, focusing on collisional [12][13][14][15][16][17][18][19][20][21][22][23] and turbulent [24][25][26][27][28][29][30][31][32] transports, although none to date was comprehensive enough to delineate nonlinear and dynamical processes occurring on actual 3D magnetic topology. An RMP can tear magnetic surfaces into resonant islands with near-resonant peeling or kinking surfaces, the latter being contributed by toroidal geometric coupling of poloidal harmonics. Previous studies have focused on the resonant response as electrons respond promptly to the island along the field lines inducing ion polarization which in turn can generate convective transport [33][34][35][36]. Isolating this polarization effect has been a significant breakthrough in capturing the RMP-induced pump-outs as extensively validated by the TM1 code [21,37,38], nonlinear two-fluid MHD modeling but on cylindrical geometry. In addition, those charged particles can also drift directly across surfaces as the toroidal symmetry is no longer retained. A numerical study [23] has been conducted on such a radial particle transport induced by Neoclassical Toroidal Viscosity (NTV) using a linear approach and found that it can considerably contribute to the pump-out. However, very few or no approaches consider these effects with nonlinear 3D magnetic topology simultaneously.
In this paper, we report an attempt to integrate these fundamental physics mechanisms into a nonlinear hybrid kinetic MHD simulation on full geometry and its validation in KSTAR experiments by reproducing the two-stage density pump-out.
As will be explained in detail, each pump-out follows distinct processes, indicating the importance of near-resonant components and integrated predictive modeling to secure access to ELM-free H-mode 3D tokamak scenarios. First, descriptions of the model and benchmark are given in section 2. Then, sections 3 and 4 present modeling results for the pump-out in two distinguished RMP-ELM suppression discharges. Lastly, the conclusion is drawn in section 5.

Nonlinear MHD response modeling
The hybrid kinetic-MHD simulations are conducted based on the JOREK code [39,40], which describes the 3D nonlinear MHD response, including the Scrape-Off Layer (SOL). Two fluid reduced-MHD equations in [41] are used for the modeling, assuming ions (i) and electrons (e) to have the same temperature (T = T i = T e ) and density (n i = n e ). Previous JOREK simulations used the ion density equation. In this modeling, however, the hybrid density evolution equation (1) is used instead of the original density equation to simultaneously include the polarization effect [42] and NTV transport.
Here n e is the electron density, e is the electron charge, S e is the electron particle source, ⃗ v e* is the electron diamagnetic, and ⃗ v ∥ is the parallel ion flow. These variables are already in the reference equations, enabling using equation (1), in JOREK. The anomalous cross-field diffusion due to microscale turbulence is represented with the empirical coefficient D assuming no changes by RMPs, similar to other Magnetohydrodynamic (MHD) simulations. Here, D has a reduced value at the pedestal to reflect the transport barrier. The first term on RHS, (∇ · ⃗ j ∥ ), represents the polarization effect in two fluids. The second term on the RHS, Γ NTV , represents particles drifting directly across magnetic surfaces due to toroidal asymmetry introduced by RMPs. This is often led by the near-resonant response and depicted as the NTV effect emphasizing its strong influence on toroidal momentum, but it can also become prominent in particle transport as suggested by recent quasi-linear single-fluid MHD modeling [23].
The simulation adopted kinetically reconstructed equilibria [43], profiles, source, radial diffusion coefficients, and plasma parameters taken before the RMP is applied from this target plasma. With the absence of the turbulence transport model, including the H-mode transport barrier, the radial diffusion coefficients and source profiles are prescribed using source balance analysis with the ASTRA-FRANTIC framework [44,45]. The coefficients on the pedestal show the values between 0.1-1 m 2 s −1 for particle, heat, and momentum channels. These diffusion and source profiles remain the same throughout the simulations. In addition, the neoclassical resistivity [46], viscosity [47], and Braginskii-like parallel heat diffusivity are employed, where the ratio between parallel and radial heat diffusivity ∼10 8 in the pedestal region. In this way, the initial plasma profile from the reference time slice can be reproduced and maintained in JOREK before applying RMPs. This approach can mimic the anomalous turbulent transport in the core and pedestal region but cannot capture its variation in time, which is a limitation of the heuristic transport model used here.
The approach in [48] is used to model the nonlinear MHD response by RMP application. We note that only n tor = 0, 1, 2 harmonics are included during the simulation to focus on the effects of RMP-induced static plasma displacement and island structures. In addition, the reference plasmas used in this study are stable for n tor > 5 edge MHD modes, enabling static RMP response simulation with n tor < 3. However, changes in edge fluctuations and ELM physics by RMP also play an important role in the pedestal transports, higher n tor > 2 components will be needed to be included in future work for more realistic simulation as they are vital for modeling edge fluctuations and ELM dynamics.

NTV modeling
NTV simulation, along with nonlinear MHD response to RMPs, has never been done prior to this work despite its potential importance since it requires information on kinetic particle trajectories ⃗ ξ deviated from toroidal symmetry. ⃗ ξ becomes the asymmetric displacement of magnetic field lines to the leading order but can be nonintegrable nonlinearly due to magnetic islands or stochastic field lines. Nevertheless, it is feasible [49] to define its radial component as ξ ψ = δT e /∇T e0 with the electron temperature in equilibrium (T e0 ) and radial temperature perturbation (δT e ), assuming the electrons along the field lines quickly reach thermal equilibrium. The figure 1 shows the displaced field lines (color contours) obtained directly from the field line tracing and the ones (blue lines) calculated using (ψ(θ) = ψ 0 + ξ ψ (ψ, θ)) with ξ ψ from the perturbed electron temperature. The reference plasma described in section 3 with I RMP =2 kA is used for these calculations. Here, ψ N is the normalized poloidal flux, and θ geo is defined as where (R 0 ,Z 0 ) is the location of the magnetic axis in the (R, Z) coordinate. As shown in the figure, this scheme has limitations in describing the detailed field structure in the island region but still shows reasonable agreement in capturing the macroscopic change of perturbed field structure and kink response around the magnetic island instead of the field structure inside the island.
Extracting ξ α , however, is limited in this nonlinear modeling as there is no effective variable showing in-surface displacement. To overcome such limitation, we approximated ξ α by linearized toroidal force balance δ ⃗ F[ξ ψ , ξ α ] = 0 based on equation (60) of [50], assuming the plasma evolves through a sequence of near-Maxwellian perturbed equilibria. Based on this approach, the in-surface displacement of field lines, ξ α , can be calculated with derived ξ ψ from temperature perturbation, which is implemented in the PENTRC code [51]. The practice deployed here to obtain both (ξ ψ , ξ α ) enabled the integration of NTV-one of the fundamental kinetic effects on toroidally asymmetric configuration, into nonlinear MHD modeling. NTV particle fluxes (Γ NTV ) are calculated from the PENTRC code, which adopts a semi-analytic drift-kinetic formulation, based on (ξ ψ , ξ α ) as prescribed. Note that Γ NTV is dominated by ions but also pumps electrons as a consequence of quasi-neutrality.
It is noteworthy that the NTV model in this work is not applicable to the inner region of magnetic islands. Nonetheless, it does capture NTV in the outer region of the islands, which can be substantially stronger than the conventional NTV without the islands. This is because, in principle, NTV scales with the strength of perturbation linearly, δB/B [52], in the neighborhood of the islands rather than quadratically, (δB/B) 2 . The importance of the NTV inside the islands must be better understood in the future, but our assumption at this point is that it could be subdominant over the other strong transports, including the parallel thermal transport, ion polarization, and E×B convection across the island flux surfaces as shown by previous work [21]. These transports alone can substantially reduce the density and temperature gradients within the island flux surfaces, so there is little room for the NTV to come into effect. This paper does not quantify nor prove this assumption, but the successful validation implies that the model might have indeed included the dominant transport processes. Furthermore, the conventional NTV formulation can be extended as long as the displacement of orbit trajectories relative to unperturbed flux surfaces can be quantified within magnetic islands. With nonlinear MHD background solutions, the displacements are not equal to the displacement of field lines nor what is called plasma displacements (which are not clearly defined unless it is ideal) but would be represented better with the equi-temperature lines based on strong parallel heat transport, which is adopted in this work. In addition, if the equi-temperature lines actually match the island structure, our representation of the lines will fail to capture all the details but still remain as a proxy to some extent, depending on the resolution of nonlinear MHD solutions.
The JOREK-PENTRC coupled simulation has been carried out in the following manner.
Step 1, impose the initial amount of RMP amplitude (Ex. 0.5 kA) in the JOREK simulation and run for 1000 Alfven time steps.
Step 2, calculate Γ NTV using PENTRC and update the Γ NTV profile in the JOREK simulation for every 1000 steps. Step 3, repeat Step 2 until the density profile and Γ NTV saturate in the iterations. Finally, in step 4, increase the RMP amplitude by 0.5-1 kA in the JOREK simulation and repeat steps 1-3. Such a discontinuous approach is to overcome the limited time scale (<1 s) of nonlinear MHD modeling. Although this is inconsistent with the experiment where RMP amplitude continuously increases over seconds, it is possible to capture the quantitative behavior of the density pedestal with gradually increasing RMP. Still, the simulation with a realistic time scale will be important for future work.

Two step pump-out modeling
This integrated simulation scheme has been used to understand the density pump-out phenomena, investigating an experimental case where the n = 1 RMP current, I RMP (or amplitude), gradually increases from 0 to 6 kA in 5 s. The studied discharge is a typical plasma reference for KSTAR n = 1 RMP ELM suppression, with the toroidal field B T = 1.8 T, plasma current I p = 0.52 MA, major radius R = 1.78 m, aspect ratio A = 3.8, elongation κ = 1.72, upper triangularity δ up = 0.38, lower triangularity δ low = 0.87, edge safety factor q 95 = 5.1, neutral beam power P NBI ∼ 3.1 MW, neutral beam torque τ NBI ∼3.9 Nm. Before the RMP application, the plasma also had the pedestal density n ped ∼ 2.25 × 10 19 m −3 , pedestal temperature T ped ∼ 0.85 keV, and normalized beta β N = 2.05. As shown in figure 2, two different pump-outs are observed at 6.5 s and 8.5 s during the increase of RMP amplitude, which corresponds to I RMP of 1.5 and 4.0 kA, respectively. Here ELM suppression is achieved at the time of the second pump-out, indicating the strong connection between pump-out and ELM suppression. One can also see the experimental pedestal density quickly recovers (in 100 ms) after the second pump-out. This can be due to changes in particle exhaustion by ELMs. Because each ELM burst expels particles to the SOL region, ELMs generate effective radially outward transport. Here, suppressing these crashes reduces particle transports, resulting in density pedestal recovery. In addition, previous work observed the change of turbulence characteristics [25] at the ELM suppression phase. Although it is unclear yet, the particle transport induced by turbulence may be changed by RMP, increasing the density pedestal.
The predicted density profile evolution resulting from increased RMP strength is presented in figure 3. Two different simulation cases are presented where Γ NTV is included or not. As introduced already, two different density pump-outs are seen at I RMP = 1.5 kA and 4.0 kA in terms of the RMP currents, as shown in figure 3(a). In figure 3(b), it is clear in the simulation that each pump-out is triggered by the rapid growth of local magnetic islands at each threshold, first by m/n = 6/1 islands at I RMP = 1.5 kA and secondly by m/n = 5/1 at I RMP = 4.0 kA. Here, both simulation cases exhibit the change of dn e /dI RMP at 1.5 and 4.0 kA with increased island size, supporting the effect of island formation on each pumpouts. At the same time, however, it can be confirmed that NTV plays an important role in the pump-out modeling, given that the case with Γ NTV is more quantitatively consistent with the experiment.
To investigate the detailed density profile evolution by RMP application, the density pedestal profile is shown in figure 4(a),  where Γ NTV is included. The figure shows the degradation of density pedestal height and gradient with increasing RMP amplitude. Here, the profile comparison between I RMP =1.5 and 2.0 kA after the first pump-out shows the density profile flattening at the 6/1 surface, supporting the idea of islandtriggered particle transport. On the other hand, the density profile change during the second pump-out (after I RMP =4 kA) exhibits a broader change instead of local flattening at the 5/1 surface. Although the profile gradient at the 5/1 surface considerably decreases, its degradation occurs over a wide range of the pedestal. This indicates the difference in detailed transport of the first and second pump-out, even though they are both driven by larger island formations.
Here it should be noted that the sudden increase of local island size above a threshold has been well understood as the field penetration process [53]. This could happen multiple times as there are multiple rational surfaces, as suggested already by cylindrical TM1 simulations. However, our simulation is rendering novel dynamics on full geometry, from the island opening to particle transport required for access to ELM suppression.
As mentioned earlier, the first pump-out results from the m/n = 6/1 island opening at the pedestal foot. However, unlike the previous two-fluid simulations with cylinder geometry, polarization effects alone cannot quantitatively account for particle transport. As shown in figure 5(a), the predicted density profile (gray line) at I RMP =2 kA without Γ NTV shows a considerable difference with the experimental profile (dotted line). This is mainly due to the weakening of polarization transport by geometry and near-resonant effects. Figure 5(b) shows the polarization term (∇ · ⃗ j ∥ ) in the electron density equation, where negative values mean density reduction. Here, the figure shows the flux averaged polarization term of equation (1), formed by the resonant component (m = 6) and near-resonant component (m = 5,7), which affects the mean (n = 0) pedestal transport. This figure shows that the cylindrical-resonant ∇ · ⃗ j ∥ at the m/n = 6/1 island is reduced by the toroidal effect (red line). In addition, the near-resonant component by edge-peeling response under the RMP with toroidal geometry strongly reduces the pump-out by polarization effect, shown as a blue line in figure 5(b).
Such weakening of polarization effect has been recovered by an additional contribution from near-resonant response with toroidicity, NTV. As shown in figure 5(a), the degraded density pedestal by polarization transport (gray line) is further reduced by NTV-induced particle flux Γ NTV (red line)  at the m/n = 6/1 island, reaching the experimental profile. Figure 6(a) shows that Γ NTV near the q = 6 surface dramatically increases when the island opens at I RMP = 1.5 kA. Here, Γ NTV peaks at ψ N > 0.99 correspond to the m/n = 7/1 peeling component. Such a drastic increase of Γ NTV on the island is due to the change in E × B rotation frequency ω E , as shown in figure 6(b). When the island opens, poloidal J × B torque exerted by the RMP is rapidly enhanced, reducing the ω E rotation near the rational surface, and Γ NTV largely increases by the drift kinetic Landau resonances effect [23,54]. A Γ NTV profile that does not take into account this braking effect, shown as a black line in figure 6(a), shows how important coupling these resonant and near-resonant effects is. Note that the NTV effects in this study are distinct from previous work such as reference [23] in that it is driven by the nonlinear evolution of ω E from island penetration rather than resistive linear response. This also explains the small effect of the 5/1 surface on the first pump-out, which has comparable island width with the 6/1 surface because of the absence of favorable ω E profile and kink response. Given that this evolution comes from island penetration, the NTV effect on pump-out can only be combined through coupling with plasma response. From this result, the strong local flattening at the 6/1 surface after first pump-out is a combined effect of island and NTV transports.
To investigate the validity of derived NTV results, we compare the measured toroidal flow and calculated NTV torque in the modeling. When we define an averaged momentum diffusivity (χ ϕ ) over the pedestal region, χ ϕ can be estimated as ∼(nV ϕ,ped − nV ϕ,sep )/∆ ped τ net , where ∆ ped is the pedestal width, τ net is net torque at the pedestal region, V ϕ,ped and V ϕ,ped are toroidal rotation at the pedestal top and separatrix. Net torque in the pedestal region (τ net ) can be approximated as τ NBI − τ NTV when NBI and NTV torques govern the momentum source. Then, assuming the increased momentum transport (χ ϕ ) with RMP and nV ϕ,ped >> nV ϕ,sep gives the following conditions, leading to the equation (4), Here, (nV ϕ ) ped,noRMP means (nV ϕ ) ped before the RMP application. Figure 8 shows the RHS and LHS of equation (4) with I RMP , where toroidal rotations of the main ion (V ϕ ) and carbon (V ϕ,C6+ ) are assumed to be similar. In addition, τ NBI is calculated from NUBEAM code [55]. As shown in the figure, the above condition is well satisfied up to I RMP < 3.5 kA in figure 7, and it can be argued that derived NTV torque and particle fluxes show reasonable values. For I RMP > 4 kA, however, τ NTV deviates from equation (4). This is mainly due to the large stochasticity of the pedestal region, which is shown in figure 7(b). The use of ξ ψ ∝ δT e becomes invalid with stochastic field lines, and derived NTV fluxes show unreasonably large values in this limit. Figure 8 thus shows the importance of the weak stochasticity assumption made in this approach while validating the model where this condition is held. In order to avoid unrealistic NTV values, we use Γ NTV profiles of 3.5 kA for JOREK-PENTRC iterations with I RMP ⩾ 4 kA. As shown in figure 3(a), the second pumpout happens at I RMP ⩾ 4 kA even without Γ NTV because it is led by the island opening. Therefore, fixing the Γ NTV profile is less likely to affect the onset of the second pump-out while may make quantitative differences, suggesting the feasibility of such an approach. However, this can still affect the plasma flow dynamics and change the island physics, raising the importance of enhanced NTV modeling in future work. Note that the NTV torque is not included in current momentum modeling and should be included in future studies because of its importance in RMP response [56]. In addition, the kinetic transport in the stochastic region is omitted here, which can be considerable with the strong stochasticity. Here, the stochasticity has a smaller effect on the trapped particle (main contribution to NTV) than the passing one, and therefore, it would have less effect on the NTV transport. However, the stochastic field is still important as it will improve the quantitative agreement between the simulation and experiment. In particular, the stochastic fields can generate effective radial particle transport through parallel diffusion. It will play a considerable role after the second pump-out (or 8 s), where strong stochasticity occurs, further degrading the pedestal gradient and height. Furthermore, recent work [57] reports the importance of kinetic transport on the pump-out in strong stochasticity. Because this modeling uses a Braginskii-like particle diffusion term where the first-order parallel particle diffusion coefficient (D ∥ = 0) is zero, it cannot treat the stochastic diffusion properly. Therefore, including the stochastic transport effect will be an important future work.
In general, the proper description of non-ambipolar transport by ion-NTV will be implementing τ NTV on the momentum equation to change the electric potential and affect the balance of the continuity equation through an electric field. However, due to the limitation of the reduced MHD equation in this modeling, direct implementation of NTV torque is difficult as the momentum equations of the reduced MHD model cannot properly treat the NTV torque as it disappears during the model reduction. Here, NTV-induced electric potential (and perpendicular flow) and non-ambipolar transport can be mimicked by including Γ NTV in the electron density equation. As the perpendicular flow is the main contributor to the rotation screening, this approach can also capture the NTV effect on the plasma response via perpendicular rotation. However, even though such NTV implementation can be effective, it still has limitations in capturing the parallel rotation change by τ NTV . Since the parallel rotation also affects the RMP screening, this model omits NTV-induced rotation screening change through the parallel direction. This limitation will be resolved by developing a full-MHD model.
The second pump-out at I RMP > 4 kA in figure 3(a) also suggests a new possible mechanism for island opening by near-resonant and mode coupling effect. It is well understood that island opening could be initiated when ω E [56,58] or ω e⊥ [19,35] profiles are near zero at the corresponding rational surface. Here, ω e⊥ is the perpendicular electron frequency. However, as shown in figure 9, both ω E (−27 krad s −1 ) and ω e⊥ (−47 krad s −1 ) for I RMP = 4 kA are not zero-crossing at the q = 5/1 surface. It turns out that the coupling between the m/n = 5/1 and 10/2 mode can enlarge the 5/1 island, where 10/2 island is a secondary island [41,59] generated by the self-toroidal coupling of near-resonant n = 1 modes. The Poincare plot in figures 7(a) and (b) show the formation of 10/2 and widening of 5/1 island through the overlap of 5/1 and 10/2 modes (figure 7(b), double O-points) when the I RMP changes from 3.5 to 4.0 kA. At the same time, pedestal stochasticity increases significantly by a spatial overlap [60] between m/n = 6/1 and another secondary island (11/2) by near-resonant mode. As a result, with an enlarged 5/1 island and stochasticity, the radial particle transport by polarization effect increases over a wide range of the pedestal, resulting in a second pump-out, as presented in figure 4(a). Thus, the nearresonant component also contributes to the second pump-out by forming the secondary mode and coupling. Interestingly, the coupling between 10/2 and 5/1 may share a similar picture with a transition to the homo-clinic topology [61]. This also may explain why the second pump-out does not exhibit strong local flattening compared to the first pump-out. It is mainly due to the broader profile degradation near the rational surface without the formation of strong Γ NTV . Still, the enhanced modeling of Γ NTV above I RMP >⩾ 4 kA will be needed to clarify it, raising the importance of future work.
It is noteworthy that the considerable deviation between the experiment and simulation occurs after the second pumpout. This can be due to the absence of ELM and turbulence dynamics discussed in the previous section. Recently, gyrokinetic studies [31,62] reported that turbulence transports also contribute to the pump-out at pedestal. Furthermore, this simulation can have considerable error bars in its quantitative results due to large uncertainties in kinetic profiles, viscosity [38], and sources. In particular, accessing the ELM suppression with the second pump-out can change SOL conditions, affecting the particle fueling and density pedestal. Complete and fully self-consistent simulations will be possible in the future with these additional mechanisms, and it will be valuable to see how these must be balanced to retain the experimental result.

Pump-out modeling in double-null like plasmas
The hybrid model is introduced to another case for further demonstration. The new application compares two doublenull-like KSTAR discharges (#29 261 and #29 270) whose parameters are consistent except for the radial width (dR sep ) of the sandwich region and upper triangularity (δ up ). The definition of dR sep is presented in figure 10(a). Their dR sep and δ up are #29 261(−1.2 cm, 0.44) and #29 270(−0.6 cm, 0.56), respectively. Their other parameters are the toroidal field B T = 1.8 T, plasma current I p = 0.54 MA, major radius R = 1.78 m, aspect ratio A = 3.8, elongation κ = 1.8, lower triangularity δ low = 0.8, edge safety factor q 95 = 5.35, neutral beam power P NBI ∼ 3.0 MW, neutral beam torque τ NBI ∼3.5 Nm. Before the RMP application, the plasma also had the pedestal density n ped ∼ 2.2 × 10 19 m −3 , pedestal temperature T ped ∼ 0.8 keV, and poloidal beta β p = 1.58.
As shown in figure 10, two-stage pump-outs are observed at 4.2 s and 5.2 s during the RMP ramp in #29 261. Here ELM suppression is achieved at the time of the second pump-out, indicating the strong connection between pump-out and ELM suppression again. Discharge #29 270 shows a similar first pump-out. In contrast, second pump-out and ELM suppression are not observed, indicating more difficult ELM suppression with smaller dR sep . A previous study [63] also reported the same effect of dR sep on the ELM suppression access. This trend has been conceptually explained with weaker RMP response and transport. However, this concept has limitations because the smaller dR sep does not lead to the weaker first pump-out, and further explanation is needed as understanding the dR sep effect is crucial in the future reactor where a double null shape (dR sep = 0) is being considered. This trend is clearer in the comparison of #29 261 (dR sep = −1.2 cm) and #29 270 (dR sep = −0.6 cm), shown in figure 11(a). Both exhibit a similar first density pedestal degradation with I RMP = 2.2 kA, while the second pump-out is observed only in #29 261 at I RMP = 3.0 kA. The integrated modeling successfully captures the characteristics of these pump-outs in both cases, as illustrated by the star markers in figure 11(a). Figure 11(b) presents the density pedestal degradation at I RMP = 2.2 kA, showing its predominant occurrence near the pedestal foot (specifically, near the q = 6/1 surface) in both cases. Therefore, this simulation confirms the role of the foot island in the first pump-out again. Here, as shown  in figure 12(a), the resonant perturbed poloidal magnetic flux (δψ pol ) at q = 6/1 surface is 45% larger in #29 261 (dR sep = −1.2 cm) compared to the other case. This finding is consistent with a previous study [63], which demonstrated that the RMP-induced plasma responses are significantly weakened by the plasma shaping effect coming from a smaller dR sep . With the presence of the resonant component at the q = 6/1 surface, both cases exhibit the formation of a foot island, with #29 261 exhibiting a 25% larger island owing to a stronger plasma response. However, figure 11. indicates that despite the difference in island size, the first pump-out level is similar in both cases.
The inconsistent trends observed in island width and the first pump-out level can be attributed to the influence of the near-resonant component. As discussed in the previous section, the near-resonant components can drive the NTV particle flux while simultaneously reducing the polarization term. These near-resonant components mainly arise from the kink response, which induces plasma displacement (ξ ψ ), as shown in figure 12(b). Here, dR sep effect also reduces RMP-induced kinking, resulting in weaker NTV in #29 270 (dR sep = −0.6 cm). Figure 13(a) presents Γ NTV comparison with I RMP = 2.2 kA, showing this overall reduction in NTV particle flux by dR sep effect. At the same time, the decreased near-resonant component contributes to an increase in the polarization effect by reducing its damping. Consequently, despite a 25% smaller island size due to the dR sep effect, #29 270 exhibits a stronger particle pump-out driven by the polarization term compared to #29 261. This observation is supported by figure 13(b), which shows the mean (n = 0) polarization term ∇ · j ∥ . Here, a more negative value indicates a stronger particle pump-out, showing that the smaller dR sep leads to an enhanced particle pump-out through the polarization term. Therefore, the polarization term and NTV in the foot island region exhibit opposing behaviors with respect to dR sep , resulting in a similar net effect in terms of the first pump-out. This example further highlights that NTV and polarization effects can be described completely through their integration and importance in unveiling the near-resonant effects.  In contrast to the first pump-out, the second pump-out exhibits considerable differences between the two cases. It is only observed in #29 261 at I RMP = 3.0 kA, accompanied by the opening of an island at the top of the pedestal (near q = 5/1), once again highlighting the role of the pedestal top island in the second pump-out. As discussed earlier, a sufficiently strong resonant field response is crucial for forming a large island [21]. In figure 12(a), the resonant perturbed field δψ pol at the q = 5/1 surface shows a 70% smaller value in #29 270 due to a weaker plasma response induced by the dR sep effect. Hence, the absence of the second pump-out in #29 270 may be attributed to an insufficient resonant response to open the island at the pedestal top. This suggests that applying a higher RMP current in #29 270 would eventually lead to the opening of the pedestal top island and the second pump-out. However, it should be noted that the maximum applicable RMP current is limited by the onset of core penetration, which can result in plasma disruption. Therefore, achieving the second pump-out in plasmas with small dR sep or double-null-like configurations becomes challenging, as the required RMP current may exceed the threshold for corelocking onset. In fact, in #29 270, the ramp-up of RMP current leads to plasma disruption before the second pump-out. This observation may explain the unfavorable effect of small dR sep on ELM suppression.
Interestingly, the similar toroidal mode coupling effect observed in the previous section is also shown in this second pump-out. Figure 14 displays the Poincaré plots of #29 261 and #29 270 for both I RMP = 2.2 and 3.0 kA. Consistent with the previous findings, a significant secondary m/n = 10/2 mode appears at the q = 5/1 surface (near pedestal top) in #29 261 case, and the widening of the top island is observed as a result of the overlap between the 5/1 and 10/2 modes. Once again, the presence of the secondary island (m/n = 10/2) is attributed to kinking or the near-resonant component. This indicates that such a toroidal mode coupling effect is not limited to previous results and suggests a non-casespecific contribution of the near-resonant component to the second pump-out. The comparison of #29 261 and #29 270 cases further highlights this point. As discussed earlier, in case #29 270, there is a reduced kinking response and approximately 30% smaller island (m/n = 5/1) structures at an RMP current of 2.2 kA, mainly attributed to a weaker plasma response and resonant components. When I RMP increases to 3.0 kA, both the q = 5/1 and 6/1 islands widen, but there is no abrupt increase in the size of the m/n = 5/1 island, as observed in #29 261. Additionally, a secondary island (m/n = 10/2) also appears, but it is significantly smaller in size compared to a larger dR sep case. The stochasticity is also considerably reduced with a smaller plasma response. Therefore, the absence of a large island opening at the pedestal top can be due to a weakened responses in #29 270, including the resonant (m/n = 5/1) and near-resonant components, which supports the possible role of mode interaction in the occurrence of the second pump-out.
It is worth noting that the m/n = 7/1 component may also play a role in the pump-outs, as it has a finite resonant component at the q = 7/1 surface (ψ N = 0.993) and a large peeling displacement, as shown in figure 12. However, the strong nearresonant (or kink) term near the separatrix substantially suppresses the polarization effect from 7/1. Additionally, the NTV particle flux is localized near the separatrix, affecting only a narrow radial region. Consequently, the impact of the m/n = 7/1 component on the density pedestal is smaller than that of the 5/1 and 6/1. Nevertheless, this component can still have a considerable effect depending on the position of the rational surface.

Conclusion
In conclusion, this work successfully captures the density pump-out by RMPs. Such bifurcating properties are led by island opening, showing a good agreement with previous results while exhibiting the difference in mechanisms by nearresonance and toroidicity. Here, the near-resonant modes play an important role by forming NTV and mode coupling. This is unveiled by integrating nonlinear MHD modeling and NTV, performed here for the first time.
Previous pump-out simulations [21,38] have been very successful in predicting the onset of first and second pumpout and ELM suppression with island physics, which is welltested in multiple devices. However, there remained questions of validity because it relies on reduced cylindrical geometry. This work does confirm again that the key to pump-out simulations is predicting the formation of the island. In that the twofluid simulation using realistic and reduced geometry exhibits similar island physics, this work explains and supports the viability of reduced geometry models in preliminary prediction on RMP physics. However, this modeling demonstrates the importance of realistic geometric effects that create quantitative differences from reduced geometric models. Here, the novel physics necessary to include in toroidal geometry in order to reproduce the same experimental observations motivates more quantitative verification of existing reduced-model predictions for future reactor scenarios. In addition, although this hybrid MHD simulation was able to obtain good quantitative agreements, there are still limitations to fully explaining the experiment. This is due to the reduced NTV modeling, absence of higher harmonics, ELM, turbulence, and source dynamics in this work. These limitations should be addressed in future work using integrated or first-principle models.
NTVs, and thanks to the KSTAR team for their experimental support. This material was supported by the U S Department of Energy, under Awards DE-SC0020372 & DE-