Modeling runaway electron mitigation via passively driven plasma instabilities in a high-current tokamak

High-current tokamaks, which magnetically confine a plasma fuel in a toroidal configuration, represent a crucial and promising next step in the quest for commercial fusion energy [1]. However, future tokamaks, such as SPARC [2], ARC [3], ITER and DEMO [4, 5, 6, 7], which will have plasma currents in the range $I_p \sim8$-20 MA, are susceptible to the formation of potentially damaging, mega-Ampere-level runaway electron (RE) beams during disruptions. This is because the multiplication of a RE seed via avalanching is exponentially sensitive to the pre-disruption plasma current, placing RE generation on power-plant-scale tokamaks in a new avalanche-dominated regime, typically inaccessible on present-day devices [7, 8, 9]. RE mitigation remains a critical, yet unsolved, challenge for fusion power-plant-scale tokamaks, where RE strikes on plasma-facing components (PFCs) may generate large, localized heat fluxes, causing melting, vaporization, or brittle failure of PFCs [10, 11], as seen on various existing machines [12].

SPARC – a high-current, high-field tokamak designed to achieve a net fusion gain $Q > 1$ [13] – will serve as a testbed for RE mitigation strategies on high-current machines [2]. Baseline massive gas injection (MGI) on SPARC will use combined neon (Ne) and deuterium (D$_2$) injection to simultaneously mitigate thermal and RE loads [3, 14, 15]. Simulations with the fluid-kinetic RE modeling framework DREAM [16], however, indicate that relatively high deuterium densities ($n_D \gtrsim10^{22}$ m$^{-3}$) are necessary to reduce the RE current below 1 MA in SPARC disruptions [17]. Moreover, a complete suppression of RE currents with material injection is generally considered difficult [9].

SPARC will therefore test a new technology – the RE Mitigation Coil (REMC) – which, as illustrated in Figure 1(a-Figure 1b), is a 3D conducting structure inside the vacuum vessel [2, 18] and represents a passive mechanism for RE suppression [19, 20, 21, 22]. Currents induced in the REMC during the disruption will seed magnetic perturbations that generate MHD instabilities in the plasma volume [19, 23]. The nonlinear interaction of these REMC-driven MHD instabilities can break the initially intact magnetic flux surfaces (Figure 1c; left), producing stochastic magnetic topology (Figure 1c; center). Stochastic flux surfaces can suppress RE growth by enhancing the transport of REs out of the plasma volume [24], with successful mitigation achieved when the loss rate limits the RE generation rate [20, 23, 25]. This typically requires preventing intact flux surfaces from re-healing in the plasma (Figure 1c; right).

While REMCs have been recently deployed on low-current tokamaks such as HBT-EP ($I_p \leq30$ kA) [26] and J-TEXT ($I_p \leq400$ kA) [27], and planned for medium-current machines such as TCV ($I_p \leq1.2$ MA) [28], SPARC ($I_p = 8.7$ MA) will be the first to test the REMC in the high-current avalanche-dominated regime. Previous computational work that combines REMC-driven MHD activity (simulated using the MHD code NIMROD) [21] with RE transport (using the particle code ASCOT5 [29]) and RE generation (using DREAM [16]) showed suppression of an otherwise unmitigated $\sim$5 MA RE beam predicted without the REMC [20]. Recent NIMROD simulations, however, show significant re-healing of flux surfaces that can confine the RE population, raising concerns about the efficacy of the REMC [21, 22].

Here, we present 3D nonlinear extended MHD + RE fluid [30] simulations of RE mitigation using the REMC that complement, and improve upon, previous disruption modeling efforts for SPARC. We demonstrate that REMC-driven MHD activity significantly reduces the RE current from multi-MA levels to <150 kA despite the re-healing of flux surfaces. The final termination of the RE beam is facilitated via scrape-off, which describes the opening-up of previously closed flux surfaces and the loss of plasma area once the plasma contacts the first wall due to vertical motion [31]. These simulations represent the first computational study that self-consistently combines 3D MHD activity, REs, and the REMC in full disruption simulations of a high-current tokamak. The simulations, performed with the extended MHD code M3D-C1 [32], offer several novel elements:

1. We self-consistently model RE generation and transport coupled with REMC-driven MHD activity, in contrast to previous efforts where different codes were used to separately model MHD activity, RE transport, and generation [20, 22], or where self-consistent MHD + RE fluid simulations of SPARC disruptions excluded the REMC [15, 33].

2. We model the full mid-plane disruption in 3D, comprising the initial fast loss of stored thermal energy called the thermal quench (TQ), and the subsequent current decay stage called the current quench (CQ). Modeling the TQ in 3D enables us to incorporate changes in the plasma current (including an $I_p$-spike), safety factor, and RE population due to MHD activity, which is important for subsequent CQ evolution.

3. We include high-Z impurity and radiative corrections to RE avalanching, providing more realistic predictions of the dominant avalanching rate, which can be significantly higher for impurity-laden disruptions [34]. Furthermore, we include fusion-relevant primary sources such as tritium beta decay and Compton scattering [33], impurity-corrected Dreicer generation [35], and prescribe a pessimistic, constant primary RE generation rate ($\sim$2 kA/ms) during the TQ to mimic the hot-tail seed reported in previous work [20, 36].

4. We use a multi-region vacuum vessel geometry, comprising the first wall, inner and outer vessels with ports, and an outer vacuum region [14, 15] (Figure 1d), as opposed to simplified vessel geometries used in previous modeling [20, 21, 33]. This multi-region geometry with realistic material-dependent resistivities allows us to self-consistently capture several important effects, such as resistive wall corrections to MHD instability growth rates and amplitudes, vertical motion, external magnetic flux diffusion, as well as currents in the vessel and the open-field-line ‘halo’ region.

5. Lastly, we model the vertical displacement and final termination of the plasma and RE beam in 3D, as opposed to previous simulations, which either did not capture this important effect in elongated plasmas [20, 21, 22, 33] or modeled the vertical displacement event (VDE) in 2D [15].

The initial equilibrium comprises the SPARC primary reference discharge, which describes a double-null diverted plasma with plasma current $I_p = 8.7$ MA and fusion power $P_{\mathrm{fus}} = 140$ MW [37]. Figure 1d shows the initial magnetic flux distribution for this equilibrium. We initialize the simulations with a pre-distributed population of roughly $3 \times10^{21}$ Ne atoms, which radiatively cools the plasma during the TQ. Although baseline MGI on SPARC will comprise combined Ne + $D_2$ injection [2], the Ne-only case simulated here represents a pathological case where we expect fast current quench rates and large multi-MA RE plateaus, as seen in previous DREAM [17],

DREAM + NIMROD [20], and M3D-C1 [15] simulations.

Figure 1b shows the simplified geometry of the SPARC REMC, which comprises two toroidal sections (major radius $R_{\mathrm{coil}} \approx2.3$ m, vertical position $Z_{\mathrm{coil}} \approx\pm0.6$ m) positioned on the low-field (outboard) side of the vacuum vessel, connected via vertical legs. In this configuration, the REMC produces a dominant $n = 1$ toroidal harmonic, and other odd $n$ harmonics whose magnitudes fall with increasing $n$ [20]. The 3D magnetic field for this REMC geometry is computed using a Biot-Savart model in M3D-C1, and the current response $I_{\mathrm{REMC}}(t)$ is calculated from the solution of a lumped circuit model

Here, $L$ and $R$ are the REMC’s self-inductance and resistance respectively, and $M$ represents the mutual inductance between the REMC and total toroidal current $I_{\phi}$ in the plasma and halo regions. The circuit model represents an improvement over previous simulations, which used a simple linear relationship between the plasma and REMC currents, and further assumed a maximum REMC current greater than the engineering limit (350 kA) for SPARC [20, 21, 22].

Results – Figure 2a shows the temporal evolution of the plasma current $I_p$ (blue), RE current $I_{\mathrm{RE}}$ (orange), and REMC current $I_{\mathrm{REMC}}$ (red dashed), with the RE current shown in more detail in Figure 2c. Here, $I_p$ includes both ohmic and RE currents inside the plasma boundary, but excludes any halo currents in the open-field-line region. A complete decay of the plasma current occurs over $\sim2.4$ ms, and the RE current remains small, reaching a maximum value of about 146 kA. The plasma current exhibits a large $I_p$-spike during the TQ, where $I_p$ rises to roughly 11 MA. As $I_p$ decays thereafter, the REMC current $I_{\mathrm{REMC}}$ increases, reaching a maximum value of about 350 kA. Instead of growing monotonically, the RE current $I_{\mathrm{RE}}$ exhibits periods of decay or slowed growth, which, as shown later, coincide with periods of MHD instability growth. When $I_p$ falls below $\sim5$ MA later in the CQ, the plasma exhibits downward vertical motion. This is illustrated by the gray dashed line in Figure 2a, which shows the vertical position of the magnetic axis as a function of time. The final termination of the plasma and RE beam is facilitated via scrape-off, as the plasma contacts the first wall near the lower divertor.

In contrast, when we artificially suppress the REMC current during the disruption, a significantly larger RE current is obtained. These results are indicated via semi-transparent lines in Figure 2(a,Figure 2c). The RE current grows to roughly 3.4 MA, more than 20× larger than that obtained with the REMC, before MHD activity during the VDE terminates the RE beam. These results, which are described in more detail in the following paragraphs, therefore indicate that the REMC plays a key role in limiting the RE current during the disruption.

Thermal Quench MHD and $I_p$-spike – Figure 2(d) shows the magnetic perturbation level $\delta B/B_0$ of different toroidal $n$ harmonics excited in the plasma during the disruption. (Note that during the TQ and $I_p$-spike, no currents were allowed in the REMC, which assumes a REMC switch.) A broad spectrum of MHD modes with dominant $n=1$ activity is excited during the TQ, with peaks in the mode perturbation level between $t \approx0.8$–$1.1$ ms. This coincides with the $I_p$-spike seen in Figure 2a. Figure 2d shows the evolution of the minimum safety factor $q_{\min}$ (gray markers) overlaid on the $\delta B/B_0$ evolution. In a tokamak, the safety factor $q(r) \equiv B_{\phi}r/B_pR$, which is a function of the toroidal and poloidal magnetic fields ($B_{\phi}$ and $B_p$ respectively) and the major and minor radii ($R$ and $r$ respectively), represents the stability of the plasma to MHD instabilities. As observed in Figure 2d, $q_{\min}$ drops from roughly 0.93 to about 0.75 initially, but jumps to about 2 by $t \approx1.1$ ms, which marks the end of the TQ MHD event.

MHD activity causes significant flattening of the current density profile $j_{\phi}(\sqrt{\psi_N})$. As shown in Figure 3a, the on-axis current density falls from about $j_{\phi} \approx12$ MA m$^{-2}$ to roughly 6 MA m$^{-2}$ at $t \approx1.12$ ms, while the safety factor profile $q(\sqrt{\psi_N})$ rises, as shown in Figure 3b, consistent with the observed flattening of current density. The internal inductance $L_i$ of the plasma decreases from about 1.24 $\mu$H initially to 0.62 $\mu$H at $t \approx1.1$ ms, acting as a driver for the $I_p$-spike. Figure 3(c–Figure 3e) shows Poincaré plots of the magnetic topology during the TQ. Tearing activity produces $(m,n)=(3,2)$, $(4,3)$, $(2,1)$ and $(3,1)$ magnetic islands in the outer $q>1$ region, breaking up flux surfaces via island overlap, while a $(1,1)$ internal kink appears in the inner $q<1$ region (Figure 3c).

As seen in Figure 3d, the magnetic topology becomes almost fully stochastic by $t \approx0.8$ ms. By $t \approx1.12$ ms, however, intact flux surfaces reform around the magnetic axis (Figure 3e), consistent with a decay in the magnetic perturbation level and a flattening of the safety factor profile after the MHD event. While the RE current grows linearly during the initial part of the TQ, as shown in Figure 2c, $I_{\mathrm{RE}}$ decreases rapidly from $\sim1$ kA to $\sim1$ A between $t \approx0.6$–$1.1$ ms, consistent with the onset of magnetic stochasticity due to TQ MHD activity. These results are consistent with previous NIMROD simulations, which also showed a rise in $q_{\min}>2$ accompanying the $I_p$-spike, together with strong losses of RE test particles during the TQ [21, 22].

Current Quench – After the TQ, the maximum electron temperature $T_{\mathrm{e,max}}$ falls to $\sim10$ eV, and remains at this low value for the rest of the simulation, providing an inductive $L_p/R_p$ time of about $\tau_{\mathrm{CQ}} \approx1$ ms. Figure 4a shows the plasma and halo current evolution in the REMC-on simulation. During the initial part of the CQ, $I_p$ decays at a rate of roughly $\mathrm{d}I_p \approx-6$ MA/ms. As shown in Figure 4a, the decay rate of the total toroidal current in the plasma-halo region is consistent with the expected exponential decay rate with characteristic time $\tau_{\mathrm{CQ}}$ (red dashed line).

A broad spectrum of MHD modes is obtained during the CQ in the REMC-on simulation, as shown in Figure 2d, but $n=1$ modes remain the dominant harmonic. $\delta B/B_0$ exhibits several peaks during the CQ, associated with the periodic growth and decay of MHD instabilities (Figure 2d). During the early CQ, MHD activity is accompanied by oscillations of $q_{\min}$ about 2. Figure 4(c–Figure 4h) shows the magnetic topology at various times during the CQ. Re-healed flux surfaces originally observed around the magnetic axis after TQ MHD activity (Figure 3e) are re-stochastized by successive MHD activity, and we recover fully stochastic surfaces again by $t \approx1.18$ ms (Figure 4c). Periodic formation and re-stochastization of (2, 1) magnetic islands is also observed during this time (Figure 4d). Intact surfaces initially present around the magnetic axis early in the CQ (Figure 3e) enable a rapid increase in the RE current, as shown in Figure 2c. By $t \approx1.3$ ms, $I_p$ has fallen to its pre-disruption value, and the RE current has increased to $I_{\mathrm{RE}} \approx3$ kA. The re-stochasticization of flux surfaces reduces the net RE growth rate, and also produces periods of intermittent RE losses, as seen in Figure 2c.

Between $t \approx1.5$ - $1.8$ ms, $q_{\min}$ falls rapidly, dropping slightly below unity. This triggers further MHD activity, producing peaks in $\delta B/B_0$ (Figure 2d), increased magnetic stochasticity (Figure 4e), and RE losses (Figure 2c) during this period. The evolution of the $q$-profile at different times during the CQ is shown in Figure 4b. The initially flat $q$-profile becomes increasingly steeper, bringing flux surfaces in the inner region into resonance with the dominant $n = 1$ perturbations driven by the REMC. The change in the $q$-profile is driven primarily by a peaking of the current density caused by a lower resistivity in the hotter, inner plasma region. Current re-distribution by MHD activity causes a subsequent rise in $q_{\min}$ after

$t \approx1.8$ ms, and the $q$-profile becomes increasingly flatter (Figure 4b). This flattening produces a region of sparse resonant flux surfaces, limiting overlap between modes on neighboring flux surfaces, and enabling re-healing around the magnetic axis, as shown in Figure 4(f-Figure 4g). Similar re-healing around the magnetic axis has been reported previously in REMC simulations [21, 22, 23]. Re-healing is accompanied by a strong increase in the RE current to roughly $I_{\mathrm{RE}} \approx146$ kA between $t \approx1.7$ ms - $1.95$ ms.

Vertical motion & termination – The magnetic topology during vertical motion of the plasma in the REMC-on simulation is visualized in Figure 4(f-Figure 4h). The plasma begins to move vertically downwards around $t \approx1.8$ ms and contacts the first wall near the lower divertor. This post-TQ plasma motion constitutes a “cold” VDE [38], triggered by the plasma current falling below a critical value, which for SPARC is roughly 5 MA [15, 39]. The plasma cross-sectional area shrinks during this time, at a rate of roughly $\mathrm{d}_t A_p \approx-2\,\mathrm{m}^2/\mathrm{ms}$, accompanied by an increase in the halo current, as current is lost from within the plasma boundary to the open-field-line region (Figure 4a).

During vertical motion, re-healed flux surfaces exist around the magnetic axis, while outer flux surfaces remain stochastic (Figure 4g). The region of intact surfaces, however, shrinks with time. As seen in Figure 2c, the RE current increases during the initial part of the VDE ($t \approx1.8$ - 1.95 ms), before decreasing afterwards. The minimum safety factor $q_{\min}$ continues to rise during this time, and the q-profile, as noted earlier, becomes increasingly flat near the magnetic axis (Figure 4b). In contrast, the edge safety factor $q_{95}$, which increases prior to vertical displacement, starts to decrease, as shown in Figure 4a. The change in $q_{95}$ is consistent with the expected scaling $q_{95} \propto a^{2}/I_{p}$, showing a decrease when the cross-sectional area begins to fall faster than the plasma current (that is, $\dot{q}_{95}/q_{95} \sim[2\dot{a}/a-\dot{I}_{p}/I_{p}]$ becomes negative). 3D MHD instabilities grow, exhibiting peaks in $\delta B/B_{0}$ (primarily even $n$ modes) around $t \approx2.3$ ms (Figure 2d). MHD activity also produces an acceleration of the RE loss rate, as seen in Figure 2c.

Comparison with the REMC-off case – A comparison of the results obtained with and without the REMC demonstrates important coupled dynamics of flux surface re-healing, vertical displacement, and RE generation. Figure 5a compares the total magnetic perturbation $(\delta B/B_{0})_{\mathrm{tot}}$ of odd $n$ (blue) and even $n$ (red) modes obtained in simulations with and without the REMC. In the early CQ ($t \lesssim1.4$ ms), the magnetic perturbation level and RE current evolution (Figure 2c) are identical in both cases. This indicates that during this period, when the REMC current is small ($I_{\mathrm{REMC}} \lesssim100$ kA), MHD activity grows primarily from perturbations pre-existing in the plasma from the TQ, instead of being actively seeded by the REMC.

Thereafter, the magnetic perturbation level becomes higher with the REMC than in the REMC-off case. The total magnetic perturbation remains high with the REMC primarily because of odd $n$ modes actively driven by the coil. This is observed in Figure 5a, which shows that the $\delta B/B_{0}$ of even $n$ modes decays significantly later in the CQ ($t \gtrsim1.5$ ms), while that of odd $n$ modes remains high. In contrast, the REMC-off case demonstrates a strong decay of both odd and even $n$ modes, consistent with the absence of any actively driven perturbations in this simulation. The lower perturbation level obtained without the REMC results in earlier and more significant re-healing of flux surfaces, as shown in Figure 5(b-c), which compares the magnetic topology between the two cases at $t \approx1.94$ ms. This is consistent with the earlier and significantly faster RE current growth between $t \approx1.5$ - 2.1 ms observed without the REMC in Figure 2c.

The substantial re-healing of flux surfaces without the REMC enables the RE population to grow and become the primary current carrier by $t \approx2$ ms. The higher $I_{\mathrm{RE}}$ slows down the plasma current decay rate, in turn, slowing down the vertical motion [^38, ^40]. This is illustrated in Figure 2a, which shows that the plasma current $I_{p}$ becomes higher and the vertical displacement is delayed in the simulation without the REMC, providing a feedback mechanism that enables the RE current to grow to a larger value without significant scrape-off losses. A similar modification of the VDE by RE currents was reported previously in 2D MHD simulations [15].

Finally, the RE current exhibits a rapid loss event (Figure 2c), associated with a large peak in the $\delta B/B_{0}$ (Fig-

During vertical motion, re-healed flux surfaces exist around the magnetic axis, while outer flux surfaces remain stochastic (Figure 4g). The region of intact surfaces, however, shrinks with time. As seen in Figure 2c, the RE current increases during the initial part of the VDE ($t \approx1.8$ - 1.95 ms), before decreasing afterwards. The minimum safety factor $q_{\min}$ continues to rise during this time, and the q-profile, as noted earlier, becomes increasingly flat near the magnetic axis (Figure 4b). In contrast, the edge safety factor $q_{95}$, which increases prior to vertical displacement, starts to decrease, as shown in Figure 4a. The change in $q_{95}$ is consistent with the expected scaling $q_{95} \propto a^{2}/I_{p}$, showing a decrease when the cross-sectional area begins to fall faster than the plasma current (that is, $\dot{q}_{95}/q_{95} \sim[2\dot{a}/a-\dot{I}_{p}/I_{p}]$ becomes negative). 3D MHD instabilities grow, exhibiting peaks in $\delta B/B_{0}$ (primarily even $n$ modes) around $t \approx2.3$ ms (Figure 2d). MHD activity also produces an acceleration of the RE loss rate, as seen in Figure 2c.

Comparison with the REMC-off case – A comparison of the results obtained with and without the REMC demonstrates important coupled dynamics of flux surface re-healing, vertical displacement, and RE generation. Figure 5a compares the total magnetic perturbation $(\delta B/B_{0})_{\mathrm{tot}}$ of odd $n$ (blue) and even $n$ (red) modes obtained in simulations with and without the REMC. In the early CQ ($t \lesssim1.4$ ms), the magnetic perturbation level and RE current evolution (Figure 2c) are identical in both cases. This indicates that during this period, when the REMC current is small ($I_{\mathrm{REMC}} \lesssim100$ kA), MHD activity grows primarily from perturbations pre-existing in the plasma from the TQ, instead of being actively seeded by the REMC.

Thereafter, the magnetic perturbation level becomes higher with the REMC than in the REMC-off case. The total magnetic perturbation remains high with the REMC primarily because of odd $n$ modes actively driven by the coil. This is observed in Figure 5a, which shows that the $\delta B/B_{0}$ of even $n$ modes decays significantly later in the CQ ($t \gtrsim1.5$ ms), while that of odd $n$ modes remains high. In contrast, the REMC-off case demonstrates a strong decay of both odd and even $n$ modes, consistent with the absence of any actively driven perturbations in this simulation. The lower perturbation level obtained without the REMC results in earlier and more significant re-healing of flux surfaces, as shown in Figure 5(b-c), which compares the magnetic topology between the two cases at $t \approx1.94$ ms. This is consistent with the earlier and significantly faster RE current growth between $t \approx1.5$ - 2.1 ms observed without the REMC in Figure 2c.

The substantial re-healing of flux surfaces without the REMC enables the RE population to grow and become the primary current carrier by $t \approx2$ ms. The higher $I_{\mathrm{RE}}$ slows down the plasma current decay rate, in turn, slowing down the vertical motion [38, 40]. This is illustrated in Figure 2a, which shows that the plasma current $I_{p}$ becomes higher and the vertical displacement is delayed in the simulation without the REMC, providing a feedback mechanism that enables the RE current to grow to a larger value without significant scrape-off losses. A similar modification of the VDE by RE currents was reported previously in 2D MHD simulations [15].

Finally, the RE current exhibits a rapid loss event (Figure 2c), associated with a large peak in the $\delta B/B_{0}$ (Fig- ure 5a), around $t \approx2.15$ ms in the simulation without the REMC, followed by eventual re-avalanching of the RE current. An inspection of the safety factor shows that the $q$-profile becomes non-monotonic due to off-axis peaking of the current profile [33]. This, in turn, triggers double tearing modes, primarily at the 2/1 and 3/1 resonant surfaces, causing strong MHD activity and rapid stochasticization of flux surfaces. Consequently, a RE plateau (i.e. a freezing of the $I_p$ decay due to the RE current) is not obtained; instead, a brief "shoulder" in the plasma current evolution is produced (Figure 2a). Similar qualitative behavior was also observed in early disruption experiments of diverted plasmas on the JET tokamak, and postulated to be a consequence of reduced RE confinement in diverted plasmas [41].

Discussion – Our results show that while a complete suppression of RE currents may be difficult to achieve on high-current tokamaks, the REMC may still provide a viable pathway for RE mitigation, even in pathological situations, as seen in these simulations. The re-healing of flux surfaces due to decaying magnetic perturbation magnitudes and a flattening of the $q$-profile following MHD activity enables a rapid rise of the RE current in the avalanche-dominated regime. However, by delaying the onset of re-healing and limiting the extent of intact flux surfaces, REMC-driven MHD activity is likely to limit the RE current to a level where vertical displacement of the plasma also contributes to RE losses during the disruption.

Recent work has explored the role of vertical displacement in limiting the avalanche gain in high-current tokamaks through reduced-order modeling of VDE dynamics [31, 42]. For vertical displacement to contribute to limiting the RE current, it is necessary that the RE current remains sufficiently low prior to vertical motion. In cases where an appreciable RE current develops, vertical motion can be slowed down (as shown in Figure 2a), and even completely avoided if a RE plateau maintains the plasma current above the critical value. This was observed when the simulations were repeated in 2D (without the REMC or 3D MHD instabilities); here, a large RE plateau of 7 MA, above the critical value, is obtained, and vertical motion is suppressed (see Figure 2b). A similar suppression of the "cold" VDE due to large RE plateaus was also observed in some 3D SPARC disruption simulations performed previously [15].

These simulations represent the first high-fidelity study of disruptions on a high-current tokamak that incorporates coupled feedback between RE physics, 3D MHD instabilities, and the REMC, together with various important disruption-relevant phenomena, such as TQ MHD activity, the $I_p$-spike, impurity effects, vertical motion, halo currents, and induced vessel currents. As noted earlier, while material injection – a RE mitigation strategy widely tested on present-day devices – can reduce RE generation rates, high injection densities and assimilation rates, which may be challenging to achieve, are necessary [3], making the REMC a valuable mechanism for mitigating REs in avalanche-dominated high-current machines. In addition to providing predictive capabilities for RE mitigation with the REMC, this work also provides valuable data for making predictions of energy deposition on PFCs, heat loads, and electromagnetic forces that inform engineering requirements for disruption-robust, power-plant-scale tokamaks.

Methods

Simulation Setup. The initial equilibrium comprises the primary reference discharge (PRD) in SPARC, which describes a double-null diverted plasma with plasma current $I_p = 8.7$ MA, toroidal magnetic field $B_0 = 12.2$ T, and volume-averaged ion density and electron temperature of $\langle n_i \rangle\approx3 \times10^{20}\,\mathrm{m}^{-3}$ and $\langle T_e \rangle\approx7$ keV respectively [37]. We note that the SPARC PRD equilibrium exhibits a $q = 1$ flux surface at roughly half the minor radius. While the SPARC PRD comprises a deuterium-tritium plasma, we model the main ion species as purely deuterium in our simulations for simplicity.

The simulation geometry comprises the first wall, inner vessel, outer vessel, and an outer vacuum region, with realistic resistivity values specified in each region ($\eta\approx8 \times10^{-7}\,\Omega\,\mathrm{m}$ in the inner/outer vessel, non-uniform parallel Spitzer resistivity in the plasma). The boundary of the outer vacuum region marks the boundary of the simulation domain and is treated as a perfect conductor. This conducting boundary is placed far away from the vessel boundaries, so as to not significantly affect the dynamics. Ports in the inner and outer vessel walls, which provide diagnostic access, are also included, and are modeled as regions of high resistivity in the toroidal direction ($\eta_\varphi\approx10^{-3}\,\Omega\,\mathrm{m}$), and low resistivity in the poloidal direction ($\eta_{RZ} \approx8 \times10^{-7}\,\Omega\,\mathrm{m}$). The domain is discretized with an unstructured mesh comprising triangular finite elements in the $(R, Z)$ plane and prism-shaped elements in the toroidal direction. M3D-C1 solves the MHD equations in the region within the first wall using reduced quintic $C^1$ shape functions in the $(R,Z)$ plane and cubic polynomial shape functions in the toroidal direction [32]. Details about the finite element treatment and numerical schemes used in M3D-C1 are provided in Ref. 32 and Ref. 43.

The Ne ionization level and emitted radiation (which comprises line, bremsstrahlung, and recombination radiation) are calculated using the KPRAD module in M3D-C1 [44]. Once ionized, Ne is assumed to be at the same temperature $T_i$ as the primary ion species. Energy is removed from the electron temperature to account for the energy required to bring the ionized Ne up to the ion temperature. Here, the injected Ne is pre-distributed in the plasma; preliminary simulations with Ne injected at MGI valve locations performed with a simplified geometry showed that the injected Ne becomes well-distributed in the plasma during the TQ. Finally, we note that baseline MGI on SPARC will comprise Ne + D$_2$ injection; this is expected to slow down the CQ time (thus, slowing down the REMC ramp-up rate) as well as the RE generation rate, as seen in previous 2D simulations which showed a maximum RE current of about 2 MA for the SPARC baseline MGI scenario [15].

Both the TQ and CQ are modeled in 3D. The TQ is modeled with 4 toroidal planes (resolving $n < 4$ toroidal modes), while the CQ is modeled with 8 toroidal planes (resolving $n < 8$ modes). The chosen toroidal resolution is sufficient to resolve the dominant $n$ harmonics excited in the system. Preliminary simulations performed with higher toroidal resolution (resolving $n < 12$) showed similar results. The primary goal of modeling the TQ in 3D is to enable MHD instabilities to redistribute the thermal pressure, current density, and impurity density, providing more practical profiles at the start of the CQ. This represents an improvement over previous modeling efforts [15, 20], where the TQ was modeled in 2D, and changes to the current profile and safety factor due to 3D MHD instabilities were therefore not included. To avoid modeling a long pre-TQ phase (which adds computational cost), the isotropic thermal diffusivity $\chi$ was initially increased to about $\sim10^3\,\mathrm{m}^2\,\mathrm{s}^{-1}$ during the TQ phase of the disruption to bring the temperature down to several-hundred eV. $\chi$ was then reduced to $\chi\sim10\,\mathrm{m}^2\,\mathrm{s}^{-1}$ by the end of the TQ.

RE Fluid Model. The RE fluid model [30, 33] in M3D-C1 is used to self-consistently calculate RE generation and transport during the disruption. Details of the RE fluid model can be found in Ref. 30. Secondary avalanche generation is calculated using the modified Rosenbluth-Putvinski formula proposed by Hesslow et al., which includes the effect of high-$Z$ impurities as well as synchrotron and bremsstrahlung radiation on the effective critical electric field and the avalanche rate, as obtained in kinetic simulations [34]. Primary RE generation via the hot-tail mechanism is expected to be significant in next-generation tokamaks [36]. The hot-tail seed is very sensitive to the TQ physics, and a direct implementation of hot-tail generation, which often requires a kinetic treatment, is not presently available in M3D-C1. Instead, we prescribe a constant primary generation rate during the TQ, providing a seed of roughly $I_{\mathrm{RE}} = 2$ kA, to mimic the hot-tail seed reported by Ref. 20 in previous simulations. We note that this value likely represents an upper bound on the hot-tail seed, consistent with predictions made using analytical models assuming fast TQ times $\sim0.1$ ms [36]. Activated sources — Compton scattering (computed using OpenMC simulations of the photon flux inside the SPARC vacuum vessel from activated walls) and tritium beta decay (assuming half the primary ion density is composed of tritium) — are also included in the present simulations [33]. Dreicer generation (both classical [45] and impurity-corrected [35]) is also available in M3D-C1 [15, 33]; here we use the impurity-corrected Dreicer model implemented as a neural network to calculate the Dreicer generation rate.

The RE fluid model used in our simulations provides the capability to tractably incorporate RE generation and transport in MHD simulations; however, a fluid treatment excludes information about electron energy, pitch angle with respect to the local magnetic field vector, and particle gyro-motion. RE losses can be sensitive to these quantities; for instance, high-energy electrons may be preferentially lost to the wall due to large drifts and gyro-orbits. MHD simulations have explored kinetic coupling [46, 47]; however, a direct coupling of kinetic REs with extended MHD codes for full disruption simulations remains challenging. Differences in RE transport and losses between the fluid model and those calculated by particle codes initialized with magnetic fields exported from M3D-C1 will be the subject of a future investigation.

REMC Model. The time-resolved REMC current is calculated by solving the lumped circuit equation with self-inductance $L = 19\,\mu\mathrm{H}$, resistance $R = 1.29\,\mathrm{m}\Omega$ and effective mutual inductance $M = 0.94\,\mu\mathrm{H}$. These values are obtained from 3D electromagnetic modeling of the REMC using the code ThinCurr [18]. The self-inductance and resistance of the coil are geometry- and material-dependent constants; the mutual inductance can change during the disruption, especially during vertical motion when the plasma magnetic axis and current centroid are displaced. A simple analytical estimate of the change in the mutual inductance $M$ between the REMC and the plasma due to changes in the current centroid $Z_c$ can be obtained by using thin coil approximations for the REMC and plasma components, and considering the fluxes through the upper and lower regions enclosed by the REMC contour (see Figure 1b) [48]. For $\lvert Z_c\rvert\le1\,\mathrm{m}$, the mutual inductance is expected to vary by $< 10\%$, representing a small change in the inductive coupling. We further note that while the plasma magnetic axis is displaced by $Z_c \approx1\,\mathrm{m}$ during the VDE in our simulation, the change in the current centroid is expected to be lower due to halo current contributions.

A Biot-Savart model calculates the magnetic field produced by the REMC in M3D-C1, and includes contributions from both the toroidal sections and the legs of the coil. Figure 6(a-Figure 6d) show the vertical and radial components of the REMC magnetic field at $\varphi= \pm\pi/2$ for $I_{\mathrm{REMC}} = 50\,\mathrm{kA}$, which are the dominant contributions produced by the toroidal sections of the coil. The field computed by M3D-C1 is in good agreement with that calculated by an independent 3D Biot-Savart solver (within 5%) and by ThinCurr (within 7%). This is shown in Figure 6(e-Figure 6g), which shows the toroidal variation of $B_R$, $B_\varphi$, and $B_Z$ at the original magnetic axis position.

The REMC’s 3D non-axisymmetric magnetic field $\mathbf{B}$ couples with the MHD equations via the scalar variables $\psi$ [where $(B_R, 0, B_Z) = \nabla\psi\times\nabla\varphi$] and $I \equiv R B_\varphi$ [43], and via the time-varying terms $\mathrm{d}_t\psi_{\mathrm{REMC}}$, $\mathrm{d}_t I_{\mathrm{REMC}} \propto\mathrm{d}_t I_{\mathrm{REMC}}$, which drive induced currents in the plasma and surrounding conducting structures. The Biot-Savart magnetic field calculation is used to construct the flux $\psi_{\mathrm{REMC}}$ and axial field terms $I_{\mathrm{REMC}}$, which are evolved in M3D-C1 instead of the magnetic field components and preserve the divergence-free condition $\nabla\cdot\mathbf{B} = 0$ [43]. The scalar quantities couple to the MHD equations, as described in Ref. 43. In addition, the time-derivatives $\mathrm{d}_t\psi_{\mathrm{REMC}}$ and $\mathrm{d}_t I_{\mathrm{REMC}}$ are added to the flux and axial field evolution equations, and induce an electromotive force in the plasma and surrounding conducting structures that screen the REMC’s time-varying field. As a simplified example demonstrating this effect, we show a simulation in Figure 6(h-Figure 6i) of a REMC-wall system, where the REMC current is first artificially ramped up to $I_{\mathrm{REMC}} \approx60\,\mathrm{kA}$, then held constant, and finally ramped down (shown by the blue curve in Figure 6h). Note that this example used a single contiguous wall for simplicity, unlike the double-walled vessel used in the simulations described in this manuscript. During the ramp-up phase, negative currents are induced in the wall, as expected from Lenz’s law. Figure 6i shows the current density in the wall during the ramp-up phase; strong negative currents are induced in the wall close to the REMC position. When the REMC current is held constant, the wall current decays resistively, and finally, when the REMC current is ramped down, positive currents are induced in the wall (see Figure 6h). A similar response is produced in the plasma region when induced currents are allowed in the plasma region (not shown here).

Computational Cost. The simulations presented in this work were performed on the NERSC Perlmutter CPU partition, and took roughly 10,000 node hours (128 cpus per node).

Data Availability

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

Code Availability

M3D-C1 is open-source software, publicly available at https://github.com/PrincetonUniversity/M3DC1.git.

Competing Interests

The authors declare no competing interests.

Acknowledgments

This work is supported by Commonwealth Fusion Systems (RPP021: Disruptions). The simulations presented in this paper used resources of the National Energy Research Scientific Computing Center (NERSC) [10.13039/100017223], a Department of Energy Office of Science User Facility (m3195 and m4829 for year 2025-26). This work is also supported by the U.S. Department of Energy under contract numbers DE-AC02-09CH11466 (MiRACL) and DE-SC0026438 (CTC). The authors acknowledge J. Chen (PPPL), C. Liu (Peking University), and E. S. Seol (Rensselaer Polytechnic Institute) for their continued support. The authors also thank J. Boguski (CFS) for help with the REMC lumped circuit model.