Toroidal modeling of runaway electron loss due to 3-D fields in DIII-D and COMPASS

In tokamak plasmas, in particular during the disruption of discharges, thermal electrons can be accelerated to a high speed close to the speed of light such that relativistic effects become important. These are called runaway electrons (REs). The acceleration is due to the presence of a large toroidal electric field generated during the thermal quench (TQ) of the plasma, when the plasma temperature quickly decreases, resulting in a large plasma resistivity, which drives the large electric field. An important factor here is that, right after the thermal quench, the plasma (Ohmic) current is still as large as the pre-disruption plasma current (in fact often even larger than the pre-disruption current during a short period of time due to a current peaking phenomenon). Consequently, a large electric field can be generated due to the Ohmic effect. More importantly, the aforementioned primary generation of the REs is often followed by a secondary generation mechanism, where close-to-speed-of-light electrons “collide” with the thermal counterparts, via the so-called large angle collision, and push the latter into the runaway regime as well. Since this “knock-on” process is proportional to the number of available REs, the secondary generation mechanism yields an avalanche effect, with a multiplication factor being proportional to the pre-disruption plasma current. The problem is, thus, particularly severe for the future tokamak device such as ITER,^1,2 which is envisaged to have as large as 15 MA pre-disruption plasma current. In fact, a large runaway current, up to 10 MA,³ is expected in ITER if a disruption does occur and the runaway beam is not efficiently avoided. These high energy, relativistic runaway electrons can be a severe threat for the plasma facing components in ITER.^4,5 For instance, a high energy RE beam may locally damage the first wall or induce a hole in the cooling pipes. This provides a strong incentive for finding adequate RE mitigation techniques in both experiments^{6–8,10–13} and theory.^14–19 Several counter-measures have been proposed to deal with the RE problem, including impurity injection (into the runaway beam)^7,13 and/or RE de-confinement utilizing three-dimensional (3-D) perturbation fields.^9–11,25 No sufficient experimental evidence has so far been obtained, demonstrating either method can robustly work in the future large scale tokamak devices. This work reports modeling results with 3-D field perturbations, either due to plasma self-generated magneto-hydrodynamic (MHD) instabilities^20–22 or via applied resonant magnetic perturbations (RMPs).^11,23–25

Significant modeling work has previously been carried out in order to understand the role of 3-D fields in the RE de-confinement. NIMROD modeling²⁰ has shown an important dependence of the RE confinement on the machine size when comparing C-mod, DIII-D, and ITER, in the presence of non-linear MHD instabilities in post-thermal quench (TQ) plasmas. Using a simplified RE drift orbit model, the runaways were, however, predicted to be well confined in ITER plasmas due to their large size. On the other hand, the R³ (R is the plasma major radius) dependence of the RE confinement time found in Ref. 20 appears to be related to the structure of the MHD perturbations during disruption, which are largely core-localized. It is unclear whether the same scaling remains if the 3-D perturbations have different structures (e.g., more edge localized as that produced by the RMP field). In terms of the RE guiding center model, Ref. 20 assumed a large aspect ratio approximation for the drift velocities of REs. While the parallel velocity of particles is time traced together with the orbit equation, a fixed particle pitch angle was assumed, which did not accurately reflect the particle's phase space variation even in the absence of any collision or dissipation (i.e., the particle's magnetic moment is not conserved). JOREK modeling²⁶ in a JET-like plasma finds that REs, initialized before the TQ, can partially survive the magnetic field stochasticity generated during the TQ, and the survival rate depends on the electron energy range. In the following work,²⁷ JOREK was applied to model the same JET-like disruption but with inclusion of an electric field. The results show RE formation due to electron acceleration by large MHD-induced parallel electric fields.

3-D RMP perturbations have been found to modify the RE confinement in small-size tokamaks^24,28 and in ITER,²³ though the latter modeling ignored the plasma response to the RMP field. For instance, Ref. 28 showed that, with the $δB/B∼10−3$ level of RMP being applied to a TEXTOR-like configuration, REs are well confined in the plasma core region but with significant loss occurring near the plasma edge depending on the particle energy and traveling direction with respect to the plasma current. Similar observations were made when the J-TEXT tokamak was modeled,²⁴ showing rapid loss of REs outside the so-called shrinkage region by the RMP field.

The present work continues our previous RE loss modeling²² based on an RE orbit module (REORBIT), which has recently been updated to allow particle tracing up to the limiting surface of the plasma facing components in various tokamak devices, including a vacuum region for diverted plasma configurations. The RE beam loss in a DIII-D post-disruption plasma, presented in Ref. 22, is re-examined with the new code capability, with results confirming our previous study and with a set of new interesting observations to be reported in this work. Furthermore, this work reports modeling results for COMPASS, where the applied RMP field has been found to be effective in reducing the runaway current.²⁵ Simulation of the RE behavior in the presence of RMP fields, for a pre-disruption ITER plasma, also reveals two different mechanisms of RE loss.

The RE loss fraction is quantified for DIII-D and COMPASS. Fast growing, plasma self-generated low-n (n is the toroidal mode number), the MHD instability is found to be the most efficient way of de-confining a well formed RE beam. The coil proximity to the RE beam and the effective coil current scaling vs the machine size and the toroidal magnetic field are two major reasons for the finite RE loss observed in COMPASS by RMP fields and the lack thereof in DIII-D. The present study particularly focuses on the RE loss pattern on the limiting surface. The lost RE distribution along the poloidal angle of the limiting surface is computed, showing a relatively wide region of the “wetted area” in DIII-D (due to the MHD instability) and COMPASS (due to RMP). In particular, REs are found to be lost to either the high field side or the low-field-side (LFS) of the limiting surface in DIII-D and COMPASS depending on the particle energy and the perturbation level.

This paper is structured as follows: Section II briefly describes the simulation models. Section III reports the RE loss modeling for DIII-D, with 3-D perturbations either coming from the MARS-F computed MHD instability or the applied RMP field. A similar study with the RMP field is also carried out for COMPASS, with results reported in Sec. IV. Section V summarizes the work.

In this work, the 3-D fields associated with MHD instabilities or the plasma response to RMP fields are computed by the MARS-F code.²⁹ The code solves single fluid, perturbed, full MHD equations in toroidal geometry. One particular aspect that we point out here is the plasma resistivity that plays an important role in all the studies presented in this work. We consistently assume the Spitzer model for the plasma resistivity across all plasmas and devices from this work.

The RE orbit tracing is launched in the presence of the total field including the equilibrium field and the MARS-F computed 3-D fields (both magnetic and electric fields). The raw representation of the computed field perturbation, in the curve-linear coordinates defined in MARS-F, is directly followed within the REORBIT module without any further coordinate transformation (e.g., into the cylindrical coordinates). This ensures the best fidelity of the simulation in the sense that the perturbed field structures, which can have sharp local variations as obtained from the resistive MHD modeling, are faithfully represented at each point in space while tracing the particle trajectory. For a divertor plasma configuration, the tip associated with the X-point of the separatrix is slightly smoothed to avoid singularity in the flux coordinates. To treat the vacuum region, a computational mesh is constructed, which smoothly connects to that in the plasma region, by fixing the poloidal angle as that defined at the plasma boundary and scaling the radial coordinate. The vacuum mesh is not an equilibrium flux surface aligned. The REORBIT module solves the following set of drift orbit equations for each individual RE,

where $(s,χ,ϕ)$ is the curve-linear coordinate system as defined in MARS-F, with $s=ψp$ (ψ_p is the equilibrium poloidal flux normalized to 0 at the magnetic axis and 1 at the plasma boundary surface) labeling the plasma minor radius, χ a generic poloidal angle, and $ϕ$ the geometric toroidal angle. $v$ is the particle's guiding center velocity to be elaborated later on. The relativistic collision timescale is $τc=mc/eEc$⁠, with $Ec=nee3 ln Λ/(4πϵ02mc2)$⁠, where n_e is the electron density, e the charge unit, $ln Λ$ the Coulomb logarithm, and ϵ₀ the vacuum permittivity. E is the parallel electric field normalized by $−Ec$⁠. The factor τ in the velocity equations is defined as $τ=τr/τc$⁠, with $τr=6πϵ0m3c3/(e4B2)$ defining the synchrotron radiation timescale. Z is the nucleus charge number of the impurity species.

Note that in the above equations, the particle momentum $P$ is normalized via $P=γmv=mcp=mc(p||+p⊥)$⁠, with m being the electron mass at rest, c the speed of light, and $γ=1+p2$⁠. The particle pitch angle is defined as $λ=(1−μ2)/b,μ=p‖/p=σ1−λb$⁠, with σ being the sign of the parallel velocity of RE. The normalized magnetic field is $b=B/B0$⁠, with B being the sum of the equilibrium field and the perturbed field and B₀ the toroidal vacuum field at the major radius of R₀. Note that the particle pitch λ here is effectively defined via the magnetic moment. The last term from Eq. (5) represents the Wiener process W_t, with increments drawn from the Gaussian distribution with zero mean and variance equal to dt.³⁰ Conceptually, W_t represents a random walk process. We emphasize that the REORBIT module chooses to time advance the particle phase space variables in terms of the total momentum and the magnetic moment, thus rigorously ensuring conservation laws in the absence of particle acceleration, collision, and radiation.

The particle guiding center velocity $v$ for relativistic REs has been derived in Ref. 31, which can be re-arranged into the following equivalent form:

where κ is the magnetic curvature and

The first terms from the right hand side of the RE velocity [Eqs. (4) and (5)] represent the parallel electric field acceleration/deceleration. The second terms are due to small angle collision, including that with the impurity nucleus, which modifies the particle pitch angle but conserves the particle energy. The third terms, associated with the τ factor, are due to synchrotron radiation. Finally, the last term from Eq. (4) represents the Bremsstrahlung drag.²² 

Equations (1)–(5) are time-advanced with initial conditions $ψp=s02,χ=χ0,ϕ0,p=p0,λ=λ0$⁠, and $σ=σ0$ at t = 0. Here, σ₀ specifies the initial direction of the RE's parallel velocity, with $σ0=+1$ defined as the direction parallel to the equilibrium plasma current (co-current) and $σ0=−1$ in the opposite direction to the equilibrium current (countercurrent). An adaptive time advance solver LSODE,³² based on the Adams predictor-corrector method, is adopted to ensure the numerical accuracy and efficiency of RE orbit tracing. The automatic choice of the time step is guided by (approximately) minimizing the certain error estimate obtained from the previous time steps and the predictor value. A conservative error tolerance of $10−7$ for the absolute error is assumed for systematic simulations. Numerical tests have been performed showing that this error tolerance level provides good convergence of the solution, with or without 3-D perturbations. Furthermore, as reported in Ref. 22, the REORBIT module has been applied to trace magnetic field lines (in the limit of the vanishing particle pitch angle and energy), recovering the equilibrium flux surfaces in the absence of 3-D fields and well-formed magnetic islands in the presence of 3-D fields in a resistive plasma.

We list new updates in the REORBIT module since our last work.²² All these updates have been included into the modeling study reported below.

• The Wiener process, i.e., the last term of Eq. (5), has been added.

• The $B*$ and $E*$ corrections have now been implemented into the code. These corrections were ignored in Ref. 22.

• The REs are now traced up to the machine limiting surface. In the previous work,²² REs were traced up to the plasma boundary, and hence, any particles trespassing the plasma boundary were considered lost. The previous treatment was not sufficiently accurate for particles located near the plasma edge since some of these particles may remain confined despite their drift orbit being extending outside the plasma domain.

• The REORBIT module now has flexible options of launching particles within any 2-D spaces out of the 5-D initial condition spaces (⁠$ψp,χ0,ϕ0,p0,λ0$⁠), representing the particle's initial minor radius, poloidal angle, toroidal angle, energy, and pitch angle, respectively. In our previous study,²² we mainly launch particles from the 2-D space of $ψp−p0$⁠, which is reasonable to study full loss of REs due to MHD instabilities. However, in certain cases (e.g., RMP) considered in the present study, the REs are only partially lost, and the loss pattern varies depending on the initial poloidal location of the particles. Therefore, we shall mainly consider launching REs from the $ψp−χ0$ configuration space in this work, while varying the particle initial energy p₀ as the third dimension.

We model the RE loss by 3-D field perturbations in DIII-D, based on discharge 177040, where a large runaway current (∼1 MA) was generated in the experiment. Large MHD bursts, with the magnitude of the kG level for the magnetic field perturbation, were observed to completely mitigate the runaway beam without regeneration in this discharge.²¹ A fast-growing n = 1 resistive kink instability (Fig. 1) has been identified as the dominant perturbation for RE suppression in this post-disruption plasma.²² In this study, we carry out further systematic investigation on the RE loss in this discharge, taking into account all the REORBIT model improvements reported in Sec. II. In particular, we highlight that the REs are traced up to the limiting surface in this work, while the RE tracing was stopped at the plasma boundary surface in Ref. 22 (due to the old REORBIT model limitation).

We start by showing the plasma shape in Fig. 1(a) and the key equilibrium profile—the safety factor—in Fig. 1(b). The beam location, being close to the center-post, was controlled in the experiment. The low-q profile, with the edge safety factor q_a just slightly above 2, is obtained for this high current runaway beam. This is the essential condition for driving MHD perturbations that grow to large amplitude within the microsecond timescale. With this q-profile, the ideal n = 1 external kink was found to be stable,²² but the n = 1 resistive kink mode, with the Lundquist number of about 10⁴ corresponding to the post-thermal quench plasma condition in this DIII-D discharge, is strongly unstable (with the growth time of tens of micro-seconds). The MARS-F computed eigenmode structure, plotted in Figs. 1(c)–1(d), shows a mixture of large $m/n=1/1$ internal kink components and the equally large $m≥2$ external kink-peeling components.

We emphasize that the MARS-F linear eigenvalue computations only identify the mode structure but not the overall mode amplitude. In the following RE loss study, we shall, thus, scan the mode amplitude, while assuming exactly the same mode structure as shown in Fig. 1. The mode amplitude, denoted as $δBp=δBpHFS$⁠, is defined by the perturbed poloidal field amplitude measured at the high-field-side (HFS) of the torus just inside the DIII-D vacuum vessel (the resistive wall). For reference, the on-axis toroidal vacuum magnetic field is 1.52 T in this discharge.

We launch REs with the same initial phase space (pitch angle and particles energy) conditions but from different positions in the configuration space. As shown in Fig. 2, the particles are initialized with a uniform distribution in the 2-D space of the plasma minor radius ψ_p and the poloidal angle χ₀ on a 30 × 30 grid. The initial toroidal angle is fixed at $ϕ0=0$⁠. At each location, two particles are launched with opposite directions for the initial parallel velocity (⁠$σ0=±1$⁠). We emphasize that most electrons have $σ0=+1$ in a typical well-formed RE beam.¹² In total, 1800 particles are launched in each simulation. We point out that a well-established RE beam does not satisfy the fluid MHD equation (the Grad–Shafranov equation). Generally, a kinetic description is needed,³³ which also means that it is more adequate to assume that the initial particles are distributed poloidally uniformly over the contour of the toroidal canonical angular momentum, not simply along the poloidal flux as in our simulations. Our approximation may not be critical for the volume integrated RE loss rate (and is certainly not important when all particles are lost at large MHD perturbation amplitude) but will affect the local properties of the RE loss. More accurate ways of launching initial REs, e.g., via specification of the particle equilibrium distribution in both the configuration and phase spaces, remain a future work.

Figures 2–6 report modeling results for REs with the initial energy of $p0=20$⁠, corresponding to about 10 MeV. The particle's initial pitch angle is always assumed to be small (⁠$λ0=0.1$ in this study), motivated by the fact that REs in well-established runaway beams are largely passing particles. The particle pitch angle of course changes during the simulation due to various scattering effects reported in Sec. II. These changes are, however, small within tens of microsecond timescale when REs get lost in order to significantly modify the particle trajectory in the configuration space.

Figure 2 shows the RE loss pattern with increasing perturbation field amplitude. The RE loss time is indicated by the color bar in these plots. For particles that remain confined, the tracing time reaches a prescribed maximal simulation time. The latter is chosen to ensure that saturation is achieved in terms of the RE loss fraction, to be reported later on.

It is evident that, without 3-D perturbation (⁠$δBp=0$⁠), co-current REs do not experience any loss [Fig. 2(a)], while countercurrent particles experience a small fraction of losses [Fig. 2(b)] for those launched near the plasma edge from the low-field-side (LFS). This insignificant loss is purely due to the large RE orbit effect in the equilibrium field. Note that the outboard mid-plane corresponds to $χ0=0° (360°)$ in Fig. 2, while the inboard mid-plane corresponds to $χ0=180°$⁠.

In general, larger RE loss occurs at increasing perturbation amplitude as expected. Nearly full loss is obtained at $δBp=103$ G, confirming our previous simulation results²² and also agreeing with experimental observation.²¹ However, different from what has been reported in Ref. 22 where REs were launched along the outboard mid-plane, Fig. 2 shows that the RE loss pattern varies depending on the initial position along the poloidal angle. With the resistive kink perturbation structure as reported in Fig. 1, the loss is not the strongest when REs are launched from either the outboard or inboard mid-plane. This pattern is more clear at moderate perturbation amplitude, e.g., $δBp=200$ G. For simplicity, we define the “loss region” as the region where the REs are launched in the configuration space and are eventually lost to the limiting surface. The loss region certainly depends on the strength of the perturbation field.

Another observation from Fig. 2 is that, with large 3-D field perturbation (close to the 1 kG level), both co- and countercurrent REs experience a similar level of loss. This is because nearly all magnetic field lines inside the plasma are lost at the large perturbation amplitude. As mentioned before, in a well-formed RE beam, most electrons travel in one direction (⁠$σ=+1$⁠). Nevertheless, we keep tracking REs traveling in both directions in this study for the purpose of comparison.

We point out that our linear eigenfunction approach is not well valid anymore, when the perturbation amplitude becomes large (say above 200 G). At large perturbation amplitude, non-linear effects generally can become important. On the other hand, our linear approximation may still be reasonable for this specific study for two reasons. First, as reported in Ref. 22, the instability grows very fast, within 10–20 μs. The RE loss also occurs within a similar timescale, as shown in the modeling (Fig. 2) and in experiments.²¹ Well-established non-linear (MHD) processes normally require a longer timescale to proceed. Second, the n = 1 perturbation has been identified as the dominant instability even when the mode grows to the largest amplitude (∼kG level).²¹ Nevertheless, we remark that a future non-linear MHD modeling of the RE loss should be valuable for this DIII-D discharge, which is beyond the scope of the present work.

The same loss pattern as reported in Fig. 2 can also be shown in the real R–Z space (Fig. 3). Interestingly, the specific “tilting” of the confined region, occurring at intermediate perturbation field amplitudes (200 and 500 G), coincides with a similar titling of the Poincare plot reported in Ref. 22 as a result of field line tracing. This shows that the RE loss correlates well with the field line loss at large perturbation amplitude.

The lost REs hit the DIII-D limiting surface. Their toroidal and poloidal locations on the limiting surface are recorded and plotted in Fig. 4 for the end time of the simulation. Certain patterns can be observed on this 2-D plane. Generally, more “striation” of the hitting points occurs with increasing field amplitude. This favorable feature is understandable since larger perturbation leads to stronger field line stochasticity, which in turns scatters the REs more. The fact that most of the REs are lost to the high field side of the limiting surface is due to two reasons. (i) The runaway beam in this DIII-D discharge is located close to the center-post of the machine and far away from the other parts of the limiting surface [Fig. 1(a)]. (ii) Most of the REs experience prompt drift orbit loss, meaning that the particles are lost to the limiting surface after a small number (sometimes even less than 1) of toroidal transits.

We emphasize that the toroidal variation of the striation pattern is obtained by launching all REs from the same (reference) toroidal angle. On the other hand, we have launched particles from different toroidal angles and found similar toroidal patterns to those shown in Fig. 4. Because the toroidal distribution of lost REs on the limiting surface is much more uniform compared to that along the poloidal angle, we shall mainly focus on quantifying the poloidal distribution of REs.

Figure 5 shows the RE loss distribution along the poloidal angle of the limiting surface after counting all lost REs along the toroidal angle but at the same poloidal location. Note that the poloidal angle is defined in terms of the arc length of the poloidal projection of the equilibrium magnetic flux surface (normalized with a periodicity of $2π$⁠), starting from the outside mid-plane. This results in the so-called equal-arc poloidal angle. The histogram is, then, translated into a distribution function by cubic-spline interpolation. Note that we normalize the peak value of the distribution to unity for the purpose of comparing different RE loss patterns at different perturbation amplitudes. The normalization factors for the distribution functions are 139, 159, and 214 for the field perturbation amplitudes of 200, 500, and 1000 G, respectively, as shown in Fig. 5(a). The corresponding numbers are 36, 107, 236, and 224 for the field perturbation amplitudes of 0, 200, 500, and 1000 G, respectively, as shown in Fig. 5(b). Generally, widening and striation of the lost RE distribution are evident with increasing $δBp$⁠, although the lost REs, with 10 MeV initial energy, still largely hit the HFS of the limiting surface. At 1 kG level perturbation (when essentially all REs are lost), the “wetted” area in the limiting surface covers about 100° along the poloidal angle near the HFS.

Next, we define the RE loss fraction by counting all particles hitting the limiting surface and dividing the total number by the number of initially launched particles. The loss fractions, as functions of the simulation time, are plotted and compared in Fig. 6, with different levels of the perturbation amplitude. We note several important effects in the RE loss shown here. First, a sharp rise of the loss fraction occurs within the initial $∼1$ μs time window, indicating prompt RE loss due to large drift orbits.²² This occurs at different perturbation levels, with stronger prompt loss at a higher field. This initial fast loss is followed by a second phase at a much slower loss rate, with the final saturation reached after about 15 μs of the simulation time. The saturated loss fraction generally also increases with the perturbation amplitude, reaching a level of 90%–100% with 1 kG field perturbation due to the n = 1 resistive kink instability.

We caution that the RE loss fraction reported in this work may slightly vary if particles were launched in different ways in the configuration space. For instance, the particles can be initially uniformly populated in the R–Z space instead of in the $ψp−χ0$ space as adopted in this work. The qualitative conclusions, however, should not be sensitive to such details. In particular, the main conclusion of nearly 100% RE loss, achieved at the 1 kG level of field perturbation with the n = 1 resistive kink structure, does not depend on how the REs are initially distributed.

We have so far assumed 10 MeV for the RE initial energy. Similar studies have also been carried out for REs with higher energy levels. Figure 7 compares the RE loss fraction for REs with initial energies of 10, 25, and 50 MeV while fixing the field perturbation level at $δBp=500$ G. The saturated loss fraction generally increases with the particle energy. The dependence is more regular for co-current REs. It is also interesting to note that, during the initial phase of simulation, the co-current REs with higher energy experience less prompt orbit loss—an observation (and explanation) also made in our previous study.²² This loss trend due to the drift effect also appears to agree with the findings from Ref. 34 where it has been shown that the magnetic transport for REs scales inversely to the particle energy in the presence of micro-turbulence. The final RE loss level, due to macroscopic field perturbations, however, increases with the particle energy, as has also been demonstrated in Ref. 28.

The RE loss pattern on the DIII-D limiting surface also varies with the particle energy. Figure 8 plots the loss distribution along the poloidal angle for 50 MeV REs. The normalization factors for the distribution functions are 449, 386, and 136 for the field perturbation amplitudes of 200, 500, and 1000 G, respectively, as shown in Fig. 8(a). The corresponding numbers are 90, 166, 227, and 256 for the field perturbation amplitudes of 0, 200, 500, and 1000 G, respectively, as shown in Fig. 8(b). At a low perturbation level (e.g., $δBp=200$ G), the lost REs are clustered in the HFS region of the limiting surface, similar to that of 10 MeV REs shown in Fig. 5. The loss distribution, however, is different from a high perturbation field (e.g., $δBp=1000$ G) for high-energy co-passing REs, as shown in Fig. 8(a). A significant amount of REs are also scattered and lost to the LFS region of the limiting surface. The total wetted area covers about $200°$ along the poloidal angle for 50 MeV co-current REs lost by the 1 kG field. On the other hand, not much REs are lost to the top and bottom regions of the limiting surface. This pattern holds for all REs with different initial energies. Similar to the 10 MeV case, the counter-passing 50 MeV REs mainly lost to the HFS of the limiting surface covering about $100°$ poloidal angle, independent of the field perturbation level.

In order to obtain a global characterization of the RE loss to the limiting surface, we define a peaking factor of the loss function as the ratio of the peak value to the mean value. In view of the machine protection, it is certainly desirable to minimize the peaking factor. Figure 9(a) compares the peaking factor of the RE loss distribution vs the resistive kink field perturbation amplitude for different assumptions on the particle's initial energy. For co-current REs, the peaking factor decreases with the perturbation amplitude. No such general trend occurs for the countercurrent REs. The difference in the peaking factor, between the co- and countercurrent REs, is largely due to the fact that the co- and countercurrent particles experience different radial drifts of the orbits. The dependence on the particle energy is generally non-monotonic, but the 50 MeV countercurrent REs generally have a low peaking factor. The lowest peaking factor of about 5, among all cases considered in Fig. 9(a), is achieved for co-current 50 MeV REs in the 1 kG perturbation field.

The dependence of the saturated RE loss fraction is more monotonic as shown in Fig. 9(b). For REs traveling in either co- or countercurrent directions, the loss fraction increases with the perturbation amplitude and the particle initial energy. Importantly, with the 1 kG field perturbation level, the loss fraction reaches nearly 100% independent of the RE's energy and traveling direction. This robust behavior is firmly consistent with the experimental observation of full loss of the DIII-D high current runaway beam by MHD bursts.²¹ 

Based on the same DIII-D plasma equilibrium, we also modeled the RE beam loss due to DC RMP fields produced by the in-vessel edge localized mode (ELM) control coils (the so-called I-coils), assuming the maximal coil current available to these coils. We found nearly no effect of the applied RMP fields on the RE loss for DIII-D. An important reason is the large distance between the I-coils and the RE beam, which is located toward the center-post of the machine. As an example, Fig. 10 shows the magnitude of field perturbation inside the plasma, produced by a 4 kAt I-coil coil current. Either the vacuum RMP field [Fig. 10(a)] or that including the plasma response [Fig. 10(b)] is more than 500 times smaller than that produced by the MHD instability in an actual DIII-D discharge 177040 inside the plasma. (We remark that the field perturbation amplitude is not the best quantity to use for characterizing the RE loss, as will be illustrated in Sec. V. Nevertheless, orders of magnitude difference in the field amplitude serves as a clear indicator on the effect on RE confinement.)

Since the DC plasma response does not modify much the vacuum field for this DIII-D plasma, as shown in Fig. 10, we also probed the idea of employing the resonant field amplification effect by applying an AC RMP field. A low-frequency response does not yield amplification since the plasma is far below the Troyon no-wall beta limit.³⁵ A scanning of the coil current frequency up to the Alfvén frequency range does find excitation of the toroidal Alfven eigenmode (TAE)-like of plasma response at the resonant frequency of 0.666ω_A for this DIII-D plasma. The amplitude of the response field inside the plasma is about a factor of 4 larger than that of the DC response, assuming the same coil current amplitude. MARS-F REORBIT modeling, with this TAE-type of response, again finds nearly no RE loss for this DIII-D plasma. We also point out that driving 4 kAt in the Alfvén range of frequency is impractical since this would require very high power voltage.

The lack of effects on REs by RMP in DIII-D is in contrast to what happens in COMPASS, which is reported below.

Figure 11 shows experimental results in COMPASS,²⁵ where application of RMP fields, produced by several kAt coil current, does result in appreciable effects on the post-disruption RE beam. In fact, Fig. 11(a) shows the increasing RE beam current decay rate with a progressive increase in the RMP coil current from zero (the reference case) to 3.5 kAt in the n = 1 configuration. The coil phasing—the relative toroidal phase of the n = 1 coil current between the upper and lower rows of the coils—is fixed at $270°$ in this series of experiments. The RE loss due to 3-D fields is evidenced by the time traces of the hard x-ray (HXR) signals shown in Fig. 11(b).

We apply MARS-F to model these COMPASS experiments. Figure 12 shows the basic equilibrium data for the modeled RE beam, reconstructed from the COMPASS discharge 15752 after thermal quench. Note that, compared to DIII-D discharge 177040, the RE beam in COMPASS is much closer to the RMP coils. As a result, the overall vacuum field perturbation amplitude inside the plasma is about 6 times larger in COMPASS, by comparing Fig. 13(a) and Fig. 10(a), with a similar coil current level (3.5 kAt in COMPASS vs 4 kAt in DIII-D).

More remarkable is the perturbation field amplitude (and pattern) including the plasma response, as shown in Fig. 13(b). The overall amplitude is nearly doubled compared to the vacuum RMP field. The non-trivial field pattern inside the plasma is due to both the $m/n=1/1$ internal kink response [q₀ is close to 1 as shown in Fig. 12(b)] and an even larger $m/n=2/1$ tearing-like resistive response, shown in Figs. 13(c)–13(d). A strong kink-peeling component in the response is also evident near the plasma edge. Modeling shows that this perturbed 3-D field structure (including the plasma response) does lead to additional RE losses as reported below, in particular for low-energy REs.

Assuming the maximal 3.5 kAt RMP coil current in the n = 1 configuration and with $270°$ coil phasing as in COMPASS experiments, REs with different initial energies are traced in the presence of 3-D fields, with results shown in Figs. 14–17. The mean average RE energy in COMPASS, calculated following a method from Ref. 36, is about 10 MeV in this COMPASS discharge. As in DIII-D, the REs are launched uniformly in the 2-D configuration space of $ψp−χ0$⁠. Figure 14 shows the increasing RE loss region with higher energy particles. For instance, substantial losses, over a large fraction of the plasma volume, occur for 25 MeV particles. The shape of the loss region is different for REs with opposite directions for the initial parallel velocity. For example, the 5 MeV co-current REs experience more loss when launched from the HFS of the torus, while the countercurrent REs experience more loss when launched from the LFS. This clear trend indicates prompt orbit loss effects due to the equilibrium field configuration, as will also be confirmed in Fig. 17.

Figure 15 shows the strike points of lost REs on the COMPASS limiting surface. Note that co-current (countercurrent) REs tend to be predominantly lost to the LFS (HFS) of the limiting surface. The exception is the co-current low energy (1 MeV) REs, to be explained later on. Counting all lost REs along the toroidal angle, at a fixed poloidal angle along the limiting surface, we again construct the distribution function as shown in Fig. 16. We notice a similar poloidal distribution in COMPASS to that in DIII-D (Figs. 5 and 8). With co-passing lost REs, the wetted region is near the HFS of the COMPASS limiting surface due to 1 MeV particles and near the LFS due to 5, 10, and 25 MeV particles. The poloidal width of the wetted area is about $100°$⁠. With counter-passing REs, the wetted area covers about $100°$ poloidal angle near the HFS independent of the particle energy.

All the observations made above are related to the difference in the RE loss mechanism in COMPASS, where the RE beam is located close to the limiting surface. As a consequence, prompt orbit loss (due to the large drift orbit width) is dominant, in particular for higher energy REs. This is further illustrated by the time traces of the RE loss fraction shown in Fig. 17. For REs with above 5 MeV initial energy, the loss predominantly occurs during the initial phase of simulation. Note that the final loss fraction at 25 MeV is very high—about 60% for co-current REs and 90% for countercurrent REs. However, only moderate loss—about 10%—is achieved for 1 MeV REs in COMPASS according to MARS-F simulation. The lower energy RE loss is less prompt, in particular for the 1 MeV co-current REs. This is in fact also the reason for the observed exception in the poloidal distribution pattern [Figs. 15(a) and 16(a)] for these particles.

It is also interesting to compare the RMP induced RE loss with the drift orbit loss due to the equilibrium field alone. The latter is not small since the RE beam is close to the limiting surface as shown in Fig. 12(a). Figure 18 compares two cases, one without RMP and the other with RMP produced by 3.5 kAt coil current. We consider REs with the initial energy of 10 MeV here. As mentioned before, this represents the average electron energy level in this COMPASS discharge. It is evident that the initial prompt loss, within several microseconds, is due to the large radial excursion of the particles' drift orbits in the 2-D equilibrium field. In particular, a fraction of particles launched from the HFS, in the co-current direction [Fig. 18(a)], drift outward⁹ and hit the limiting surface at the LFS [Fig. 18(e)]. Conversely, a fraction of particles launched from the LFS, in the countercurrent direction [Fig. 18(c)], drift inward and hit the limiting surface at the HFS [Fig. 18(g)]. This prompt loss due to 2-D fields occurs within a single poloidal transit and quickly saturates in time. Inclusion of RMPs produces more RE losses at longer timescales (hundreds of microseconds). The loss patterns on the limiting surface remain similar to those due to equilibrium fields.

Figure 19 quantifies the loss fraction and the loss pattern. Compared to the equilibrium field induced drift orbit loss, we find that the RMP field increases the loss fraction by about 70% for the co-current REs and about 50% for the countercurrent REs in this COMPASS discharge. Note that the loss saturates after about 200 μs [Figs. 19(a) and 19(b)]. The wetted area on the limiting surface, shown in Figs. 19(c) and 19(d), mainly covers the LFS due to co-current REs and the HFS due to countercurrent REs. Application of the RMP field only slightly modifies the RE loss distribution.

Figure 20 shows two typical examples of the RE trajectories with lower and higher energies. With both particles being launched from the same location (HFS), the 5 MeV co-current RE [Fig. 20(a)] experiences fast loss to the outboard side of the COMPASS limiting surface (thick black line). This is essentially due to the large drift orbit of the RE in the equilibrium field. The 3-D perturbation field plays a minor role. On the contrary, the co-current 1 MeV [Fig. 20(b)] RE is confined without the RMP field. Nevertheless, the presence of the 3-D field modifies the particle orbit, leading to eventual loss to the HFS of the limiting surface after many transits along the poloidal angle (and the toroidal angle).

It is interesting to elaborate on possible reasons for the qualitative difference between the RMP induced RE loss in DIII-D (no loss) and COMPASS (moderate loss). One reason has been discussed before, i.e., the much further distance between the RE beam and the RMP coils in DIII-D, leading to ∼6 times weaker perturbation field in the plasma region as compared to that in COMPASS. The plasma response induced field amplification in COMPASS further enlarges this difference in the perturbation amplitude by a factor of about two. Another important reason is the machine size and the toroidal magnetic field strength. It is easy to show via dimensional analysis that what matters for the RE confinement is the effective RMP coil current $Ieffective=IRMP/(R0B0/μ0)$⁠, where $IRMP$ is the applied coil current in kAt, R₀ the plasma major radius, and B₀ the toroidal vacuum magnetic field defined at R₀. For DIII-D discharge 177040 at 1025 ms, we have $B0=2.58$ T at $R0=1.4$ m. For COMPASS discharge 15752 at 1080 ms, we have $B0=1.07$ T at $R0=0.57$ m. Therefore, even with the same coil current in kAt, the effective current in COMPASS is about 6 times larger than that in DIII-D.

We point out that synthetic diagnostic data are needed in order to enable direct comparison of the RE loss modeling results with the HXR data as shown in Fig. 11(b). Constructing such a synthetic diagnostic, as well as further developing RE launching options with given distribution functions in both the particle configuration and phase spaces, remains a future work. Finally, establishing scaling laws vs the machine size and/or the RMP perturbation strength requires more systematic modeling efforts, which is beyond the scope of the present investigation.

The 3-D field induced RE loss has been simulated for DIII-D and COMPASS plasmas, utilizing the MARS-F code incorporated with the recently developed and updated REORBIT module. The key findings are summarized in Table I.

In the DIII-D discharge 177040, the post-disruption beam is found to be nearly completely (>90%) de-confined, by 1 kG level of field perturbation, produced by a fast growing n = 1 current driven resistive kink instability occurring when $qa∼2$⁠. The loss time is about 15 μs. We note that the exact value of loss fraction here slightly varies depending on the initial conditions for REs. The particles are uniformly launched in the configuration space in this study. In our previous study,²² 100% loss was computed in a mixed space of particle minor radius and energy. These basic results are consistent with experimental observations.²¹ 

On the other hand, among several REORBIT model updates compared to that reported in Ref. 22, one important new aspect is the extension of the RE tracing up to the DIII-D limiting surface (instead of assuming the plasma boundary as the limiting surface). Furthermore, the previous study mainly considered the RE launching position along the LFS mid-plane, while this study finds non-uniform loss patterns (loss region) when particles are launched from different poloidal locations inside the plasma. In fact, the loss region is found to correlate very well with the region of the magnetic field line loss in the Poincaré space. The largest loss region is neither in the HFS nor LFS of the plasma domain.

Much work is also focused on the RE loss pattern on the limiting surface in this study. A loss distribution function is defined along the poloidal angle of the limiting surface after counting all lost REs along the toroidal angle. The loss distribution function shows the poloidal location and coverage (the wetted area) of lost REs on the limiting surface. Counter-current lost REs typically hit the HFS of the limiting surface, covering about $100°$ poloidal angle. For co-current REs (which consist of the dominant portion in a well-established RE beam), the wetted area depends on the particle energy. While low-energy REs are still lost to the HFS of the limiting surface covering about $100°$ poloidal angle, the high-energy REs can also hit the LFS covering about $70°$⁠. More poloidal striation is observed in the loss distribution function, with either the increasing perturbation field amplitude or the increasing particle energy.

Another quantitative measure for the wetted area is the peaking factor of the loss distribution function. Because the loss function is always peaked in a certain poloidal location, the peaking factor is always relatively large. For 50 MeV co-current REs in the 1 kG level perturbation field, the peaking factor can be as low as 5.

For the same RE beam in DIII-D, we find that application of the DC RMP field, produced by 4 kAt I-coil current, does not help to de-confine REs. Almost no RE loss is found in the modeling. Another idea has also been pursued with the RMP field by scanning the excitation frequency up to the Alfvén range. A TAE-type of plasma response is indeed found at the excitation frequency of 0.666ω_A for this DIII-D plasma, where the perturbation field amplitude is increased by about a factor of 4 compared to the DC response. This TAE-type of perturbation field, unfortunately, still does not produce any appreciable loss for this DIII-D RE beam.

The situation is different in COMPASS where a finite fraction of RE loss is obtained in the presence of the DC n = 1 RMP field produced by 3.5 kAt coil current. About 10% loss fraction is modeled for the low-energy (∼1 MeV) REs. The loss fraction reaches up to 90% for high-energy particles (25 MeV). On the other hand, this large loss in COMPASS is shown to be largely due to the finite drift orbit width effect, which occurs (at a slightly lower fraction) even in the absence of the 3-D field. An important reason for the large drift orbit loss is the proximity of the RE beam to the COMPASS limiting surface. Poloidal striation is also observed in the loss distribution function for higher energy REs, similar to that in DIII-D.

Several factors contribute to the qualitative difference in the RMP induced RE loss between DIII-D and COMPASS. One is the proximity of the RMP coils to the RE beam. This distance is much larger in DIII-D. As a result, with a similar level of the applied coil current, the vacuum field is about six times smaller in the plasma region in DIII-D, compared to that in COMPASS. Second, a factor of about 2 in plasma response induced field amplification is obtained in COMPASS, while no field amplification occurs in the DIII-D plasma. Finally and importantly, what matters for the RE confinement, as long as the 3-D field perturbation field is concerned, is the effective RMP coil current, which is inversely proportional to $R0B0$⁠. The smaller machine size and lower toroidal field in COMPASS contribute another factor of 6 in favor of RE de-confinement in COMPASS.

The above scaling is certainly not favorable for ITER. On the other hand, with the $∼10$ MA level of the post-disruption RE beam current as predicted for ITER (due to the extremely large avalanche amplification factor), it is expected that the low-q plasma is achievable after the thermal quench. Therefore, the very efficient RE beam mitigation scheme by fast growing MHD events, as demonstrated in DIII-D, may be one of the potential solutions for the RE problem in ITER, provided that the RE seeding and avalanche cannot be avoided after the thermal quench. The RMP field, however, may be helpful to change the RE seeding in ITER, which remains a future investigation.

As a final remark, we list important factors that affect the RE loss fraction and the wetted area in the limiting surface based on the findings from the present study: (i) the perturbation structure and amplitude, in particular a high amplitude perturbation, which is critical for achieving large and uniform (in both configuration space and particle phase space) RE loss; (ii) the particle energy, typically with higher saturated loss for higher energy REs; (iii) the proximity of plasma to the limiting surface, which (together with the perturbation structure) affects the poloidal location of the wetted area and the loss mechanism (e.g., the drift orbit loss).

This work was supported by the U.S. DoE Office of Science under Contract Nos. DE-SC0016452, DE-FC02–04ER54698, and DE-FG02–95ER54309. This report was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof. This work was supported by the Czech Science Foundation (GACR) under Grant No. GA18–02482S and the MEYS projects LM2015045 and CZ 02.1.01/0.0/0.0/16_019/0000768 and carried out within the framework of the EURO-fusion Consortium. This work was also supported by the EUROfusion—Theory and Advanced Simulation Coordination (E-TASC). It also received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement No. 633053 co-funded by MEYS Projects Nos. 8D15001 and LM2018117. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

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