Boundary condition effects on runaway electron mitigation coil modeling for the SPARC and DIII-D tokamaks

Extended-MHD modeling of planned Runaway Electron Mitigation Coils (REMC) for SPARC and DIII-D is performed with the NIMROD code. A coil has been designed for each machine, with the two differing in shape and location, but both having n = 1 symmetry (with n the toroidal mode number). Compared to previous modeling efforts, three improvements are made to the simulations boundary conditions. First a resistive wall model is used in place of an ideal wall. Second, the ThinCurr code is used to compute the time-dependent 3D fields used as magnetic boundary conditions for the simulations. Third, the simulation boundary is moved from the first-wall location to the Vacuum Vessel (VV), which extends the boundary past the location of the internal REMC. To remove the 3D coil from the simulation domain, an equivalent set of 3D fields is calculated at the VV boundary that produce approximately the same field distribution at the last closed flux surface assuming vacuum between the two. Each of these three boundary condition improvements leads to an improvement in the predicted performance of the REMC for both machines. The resistive wall alone primarily effects the resonance of the coil with the plasma after the TQ, affecting the q-profile evolution in the SPARC modeling, and allowing the applied spectrum to be modified in response to the plasma in the DIII-D modeling. The movement of the simulation boundary has the most significant effect on the RE confinement overall, including in the early stages, particularly for a DIII-D inner wall limited equilibrium, where the RE loss fraction increases from 90% to > 99%, with SPARC RE losses also occurring much earlier when the boundary is placed at the VV.


Introduction
The present class of experimental tokamaks can not replicate all of the disruption consequences for a reactor-class tokamak, nor can they provide a sufficient basis to ensure disruption-free operation over the multi-decadal life-span of a future reactor.A tokamak fusion power plant can not be constructed under the assumption that no disruptions shall occur; rather, disruption resiliency must be incorporated into the design of a tokamak reactor [1].Active disruption mitigation strategies based on massive material injection will likely be incorporated into a disruption resilient tokamak power plant design, but passive solutions are especially attractive.The Runaway Electron Mitigation Coil (REMC) [2][3][4][5] is one proposed passive solution to the problem of disruption generated runaway electrons (REs) [6].The concept relies on the large loop voltage induced during the current cuench (CQ) phase of a disruption to drive current in a 3D coil in order to produce large 3D magnetic perturbations to promote RE orbit losses faster than REs can be generated.
The dangers to reactor integrity associated with REs of 10s of MeV energy carrying multiple MAs of current have not been resolved largely because the problem has not been entirely realized in present day tokamaks.Disruptions, in which a loss of confinement leads to a thermal tuench (TQ) followed by a CQ of the highly resistive cooled plasma, occur routinely in present tokamak experiments, often without significant threat to the machine.Nevertheless, numerous instances of RE induced wall melting have been observed even in present machines [7][8][9].A reactor will differ from tokamak experiments in two primary respects.First, the RE 'knockon' avalanche [10], which causes exponential growth of a seed RE population due to small-angle collisions between REs and thermal electrons, has a total gain that is proportional to the exponential of the initial plasma current-placing ∼10 MA tokamaks in a qualitatively different regime of avalanche multiplication than ∼2 MA tokamaks.Second, nuclear devices will have seed sources from Compton scattering and tritium beta decay [11] that are continuous and ineliminable.The combination of these two factors means that attempts to control or deconfine the hot-tail RE seed population produced during the TQ are effectively irrelevant to the final RE current [12], so that strategies like two-stage SPI injection in the TQ are at best stop-gap measures that may ensure success of the nonnuclear phase of ITER but do not translate to a reactor.In a reactor, the problem occurs in the CQ and the solution must be in the CQ.Absent a strategy to induce sufficient RE losses in the CQ, only the 'benign termination' strategy of low-Z injection into an already formed RE plateau has shown promise for preventing catastrophic machine damage [13,14].These important recent results on JET and DIII-D represent the only demonstrated strategy applicable to ITER nuclear operation, but ideally this would be a plan of last resort in the case that RE plateau prevention fails.
The REMC concept leverages the problem (the high CQ loop voltage) into a solution (current driven in a 3D coil) so that it does not depend on the success of a disruption prediction algorithm to achieve RE mitigation.The concept has not been experimentally tested to date, but coils have been designed for both SPARC [4] and DIII-D [5] and should be operational in one or both devices by 2026.While DIII-D has an avalanche gain factor of only 50-100, it is capable of producing postdisruption RE currents of several hundred kAs and is a flexible experiment [15] with excellent diagnostics for REs [16].The SPARC device [17], which is still under construction, will be a mission driven experiment but has a reactor relevant avalanche gain factor of around 3 billion and will operate with deuterium and tritium.The two coils have different designs, but each has predominantly n = 1 symmetry (with higher-n modes also driven at lower amplitude), where n is the toroidal mode number.The DIII-D coil has a single helical turn along the inboard wall connected by a single vertical leg [5].The SPARC coil has upper and lower toroidal segments on opposite sides of the outer wall, connected by two vertical legs [4].The geometry of the two coils is shown in figure 1. Results from these two experiments combined will provide a firm basis for model validation and qualification of the concept for future Fusion Pilot Plant (FPP) designs.

SPARC modeling
Three designs for the SPARC coil having predominantly n = 1, 2, and 3 symmetries, respectively, were initially proposed [4].The three designs were evaluated [19] using a combination of electromagnetic response modeling with COMSOL and nonlinear magnetohydrodynamic (MHD) modeling with NIMROD [20].The maximum coil current and 3D field distribution in response to a prescribed linear current ramp-down were used to apply time dependent 3D field boundary conditions in NIMROD.The NIMROD simulations also included test-particle drift-orbit calculations for REs [21] as the MHD fields evolved.In the three simulations, calculated RE losses for the n = 2 and n = 3 coil were found to be clearly insufficient to overcome the high SPARC avalanche gain, while further modeling was pursued to fully qualify the n = 1 coil.
For a more complete evaluation of the n = 1 coil, two additional codes were included in the workflow [22].The ASCOT5 code [23] was used to determined RE transport coefficients for many individual time slices using NIMROD 3D fields.For this calculation, a larger number of drift-orbits was integrated at each time to evaluate advection and diffusion coefficients over a range of energy, pitch and minor radius.These RE transport coefficients were then used in a RE generation and evolution simulation with the DREAM code [24].With the CQ calculated independently in DREAM, the transport coefficients from ASCOT5 were mapped based on the same value of plasma current in NIMROD.The DREAM modeling without the effects of the REMC found a maximum RE plateau current of ∼5 MA, whereas with the transport from the coil included, full prevention of the RE plateau was achieved [22].
Further simulations performed for the SPARC n = 1 coil focused on the effects of including the TQ MHD in the NIMROD modeling [19], and the effects of varying assumptions about transport inside small islands in ASCOT5 [25].It was found that when the MHD activity associated with the TQ was included in the simulation, earlier and faster RE losses were obtained in this phase.However, the current profile redistribution associated with plasma current spike at the end of the TQ was also found to lead to a rapid increase in the onaxis safety factor, q, and once q exceeded 2 everywhere, a core of good flux surfaces reappeared sooner than in the simulation including only the CQ [19].The effects of confinement in small remnant or reformed core islands was explored with a set of DREAM simulations in which these assumptions were varied.Although the levels of transport obtained directly with ASCOT5 lead to full RE suppression, if instead zero transport in these islands was assumed, up to ∼1 MA of RE current was obtained.In a scan, the minimum level of transport needed in the islands to recapture the full prevention result was found to be well below the values obtained with ASCOT5 [25].Further simulations examined the effects of the close ideal wall, which was placed at the first-wall (FW) location in order to use the fields from the internal coil as boundary conditions.A comparison of TQ-only simulations (no coil) with a close and far wall showed a significant stabilizing effect of the FW boundary condition, suggesting this may also be an important effect in the REMC simulations [19].

DIII-D modeling
The DIII-D coil was placed along the inboard wall because CQ plasmas are typically inner-wall limited (IWL).The design was down-selected based on Vacuum Island Overlap Width (VIOW) calculations [5], and two candidate designs were modeled with the linear-MHD response code MARS-F [26].In these simulations, equilibria from the mid-CQ phase of a disruption were used and the linear response to the vacuum fields from an assumed half-amplitude coil of 50 kA or 100 kA was calculated.The loss fraction of REs launched uniformly in the poloidal flux coordinate was calculated for each case, and up to 70% loss was found depending on the q-profile and coil current [5].
Nonlinear NIMROD modeling for the chosen inner-wall coil design was performed using the vacuum fields as boundary conditions and assuming a linear response of the coil current to the plasma current decay, with maximum coil currents of 100 kA or 200 kA, corresponding to the mid-CQ values from the linear modeling.Beginning with an IWL early-CQ equilibrium, the nonlinear response resulted in losses of ∼ 90% of test-particles, almost independent of the maximum REMC current.The surprising insensitivity to the coil current amplitude was explained by the fact that the initial island overlap consistently occurred when the edge-q exceeded 8, triggering an inward cascade of stochasticity.Beginning instead with a lower-single null (LSN) diverted equilibrium, an increasing loss fraction with increasing coil current was found, with losses up to 98% of the initial population.

Limitations of the prior modeling
In preparation for validation of REMC simulations, improvements to simulation fidelity are being pursued, particularly with respect to boundary conditions that adequately represent the effects of an internal 3D coil located between the FW and the vacuum vessel (VV).Difficulties in this respect arise because NIMROD can not accommodate a 3D conducting structure within the gridded simulation domain, leading to the placement of the simulation boundary at the FW location.The ideally-conducting wall boundary conditions used at the FW location in the prior modeling can have a particularly stabilizing effect on MHD modes triggered by the coil (or the TQ), but substitution of resistive wall boundary conditions at the FW only partly addresses the problem, where one would really like to move the resistive wall to the location of the experimental VV.
The NIMROD code has multiple implementations of a resistive wall boundary condition that have usefulness in different physical situations.Both implementations are a thin resistive wall model, which assumes a jump in the magnetic field at the wall ∆B between the inner plasma solution and the outer vacuum solution.The jump corresponds to a surface electric field that produces a change in the normal magnetic fields at the wall.
Here the wall velocity v wall = η wall /µ 0 δ wall , where η wall is the wall resistivity and δ wall is the wall thickness, is the single parameter specified in NIMROD to characterize the wall.The jump in the magnetic field requires a solution for the external vacuum fields, and this can be obtained with a solution to vacuum field equations in a separate gridded region outside the resistive wall.This method is not used here, but its use in NIMROD is descried in [27].Another implementation with no external gridded regions makes use of a Green's function method for coupling the vacuum field at boundary points.The Green's function implementation is more straightforwardly compatible with the present placement of the boundary inside the coil, but the present version does not include the n = 0 component of the resistive wall.The implementation in the n > 0 components is simplified by the NIMROD Fourier representation in the toroidal direction, which allows analytic integration toroidally.Moreover, the n > 0 components of the normal magnetic field inherently integrate to zero over the bounding surface.For the n = 0 component, if the change in the normal magnetic field at the boundary is not obtained with enough accuracy, the total flux into or out of the volume may be non-zero, forcing a magnetic monopole in the simulation domain.These numerical issue should be surmountable, but for now the n > 0-only boundary condition is used.The consequence is that axisymmetric modes will produce an ideal wall response.For instance, growth of vertical instability is precluded with the n > 0-only resistive wall model.
Modeling for both machines in the subsequent sections is performed for the chosen n = 1-symmetry coil design for each machine.The plotted spectra of toroidal mode amplitudes imposed by the predominantly n = 1 coil or destabilized in the simulation include toroidal mode numbers n = 1 − 10 (the limit of the simulation resolution).We will first compare idealwall-at-the-FW simulations to the same cases including the n > 0-only, Green's function resistive wall model for the same geometry.We will then describe two further improvements to the modeling: the incorporation of results from the ThinCurr electromagnetic response code, and extrapolation of the fields at the last closed flux surface (LCFS) to a larger boundary shape.Individually and in combination, these model improvements lead to enhanced RE losses compared to the results in [19], indicating that preliminary modeling was in many respects pessimistic in predicting the REMC effectiveness in either device.

SPARC comparison
As stated, the resistive wall in NIMROD is specified by a single parameter v wall which has units of velocity and is equal to η/δµ 0 , where η is the wall resistivity and δ is the wall thickness.For the SPARC resistive wall simulations, we use a value of v wall = 40 m s −1 based on just the inner of two conducting walls on SPARC.The resistive wall boundary condition is applied to the n = 1 − 6 toroidal Fourier components, where the resistive wall for the n = 0 component is not available, and maintaining an ideal wall for n = 7 − 10 is found empirically to aid in numerical stability (for reasons that would require further investigation).In this section, we repeat an ideal wall simulation that was presented in [19], which included MHD associated with a TQ triggered by Ne injection followed by the typical I p -spike seen just before the CQ phase of a disruption, and the subsequent CQ.In the ideal wall simulation, the I p -spike was followed by a redistribution of the current profile and a rapid increase of q 0 .Coinciding with q 0 < 2, two islands appeared and then rapidly coalesced into a large region of closed flux in the core, which persisted through the CQ phase.
In the resistive wall simulation, a slightly earlier TQ leads to an earlier but comparably sized I p -spike, followed by (most importantly) a much gentler rise in the q 0 after the current spike (figure 2).After the TQ, the flux surfaces in the core do transiently reappear.However, once the REMC amplitude grows further in response to the current decay, these core flux surfaces are again stochasticized, unlike the ideal wall case.This situation similarly persists only until the on-axis q rises above 2, but this occurs at a much later time relative to the overall current decay.With the avalanche gain ≈exp(2.5Ip [MA]) [28], every additional MA of current that is dissipated before good confinement of REs is restored can reduce by more than an order of magnitude the RE current multiplication that can occur in the confined region.

DIII-D comparison
Whereas previous modeling for DIII-D included only the CQ phase of the disruption and the plasma response to the coil, here we compare ideal and resistive wall simulations that include a TQ induced by Ne radiation (figure 3), resulting in stochastization of the magnetic field prior to the start of the CQ.Here, we deposit Ne across the entire plasma volume at the same time, rather than trying to mimic a particular Ne injection scheme in detail.As a result, the TQ occurs in two stages, with the initial rapid cooling by ionization, radiation and dilution followed by an MHD-triggered final phase of the TQ.The global evolution of thermal energy and plasma current is not much affected by the ideal vs. resistive wall, but the TQ ends earlier when the REMC is active due to the coil With a resistive wall (lower row), the SPARC disruption has a more gentle increase in q 0 after the Ip-spike, compared to the ideal wall simulation (upper row).In both simulations, some flux surfaces reform in the core shortly after the TQ, but as the REMC grows during the CQ, the core flux surfaces are again destroyed when q 0 remains below 2.
response induced after the first phase of the TQ.A population of 10 530 RE test-particles is launched at t = 0, having a radial distribution that is uniform in normalized poloidal flux, and energy distribution that is uniform between 1 and 30 MeV, and a range of pitch p ⊥ /p that is uniform from 0-0.5.This choice does not mimic any particular expected distribution, but rather allows for a simple evaluation of preferential loss from particular regions in real or momentum space.The energy of the individual REs is not evolved in the simulation, so that the behavior of lower energy REs can be assessed even late in the simulation.In each case, the initial test population is lost during the late TQ, and a new test-population of the same size is launched 0.5 ms after the initial loss is complete.This time interval is simply chosen to allow sufficient reformation of closed flux that some of the new test particles are confined.Here we see noticeable differences in the confined RE fraction among the four cases.Both resistive wall cases have a larger fraction promptly lost than the corresponding ideal wall case, but the largest separation is between the resistive wall simulation with the REMC and the other three.Whereas the other three simulations have around 98%-99% of test particles promptly lost, the resistive wall simulation with the REMC promptly dumps closer to 99.8% of REs.In other words, the combination of the REMC and the resistive wall has a larger effect on RE confinement than either of the two effects individually.
During the CQ, we see that the growth of the n = 1 magnetic energy driven by the REMC is faster for the ideal wall case.On the other hand, the relative (to the n = 1) amplitude of the n = 2 component is larger for the resistive wall case.With the ideal wall, the field lines are strictly frozen into the boundary and the externally applied amplitude of each component is fixed.With a resistive wall, the field lines can slip along the boundary and redistribute both poloidally and toroidally in response to the plasmas, so that, for instance, energy can be transferred from the n = 1 component to higher n components, as happens here.
We can see this clearly if we compare the toroidal mode energy spectra late in time when the external fields strongly dominate (figure 4).With a resistive wall, the energy in the n = 1 mode is substantially reduced while the n = 3 − 6 components are all enhanced (with n = 2 being almost unchanged).Because the ideal wall boundary condition is maintained for the n = 7 − 10 components, those match the ideal wall simulation.
The larger loss fraction of REs with the resistive wall is correlated with a delayed reappearance of good flux surfaces in the core (figure 5).At 1.3 ms, a small well-confined region reappears near the axis only in the ideal wall simulation.At this time, neither simulation has a q = 1 surface in the plasma, and the low-order rational surface nearest the axis is 3/2.In the ideal wall case, where the applied spectrum at the boundary is fixed, the 2/2 component is much larger than the 3/2 component.As with the toroidal spectrum, the resistive wall case has a much flatter poloidal spectrum and consequently a larger 3/2 component.By 1.5 ms, both simulations have formed 2/1 islands, but the axis remains stochastic in the resistive wall simulation.
An additional effect of the boundary condition is observed in the pattern of RE test-particles striking the simulation boundary.In the resistive wall simulation, REs strike the boundary at only two poloidal locations (covering 360 • toroidally): namely, the outer divertor strike-point and the outermidplane.Lower energy REs follow the magnetic field lines to the divertor, while higher energy REs strike the outer wall due to curvature drift.In the ideal wall simulations, some REs strike other locations.A few lower energy REs strike just above the inner mid-plane, with a toroidal extent of 240 • , while a few high energy REs strike the upper divertor with a toroidal extent roughly corresponding to the gap in the innermidplane strikes.These striking locations are attributed to the localized radial fields produced by the coil along the inner wall, which can divert some REs orbiting near the LCFS into the wall.It is observed that the amplitude of these localized fields is attenuated (by around a factor of two) in the resistive wall simulations.Again this is due to the ability of the applied fields to spread out along the wall in response to the plasma.

Electromagnetic response modeling with ThinCurr
For the preceding and previously published [19,22] simulations, the distribution and ramp rate of 3D fields produced by the coil were calculated using COMSOL for SPARC, and using vacuum fields with an assumed linear ramp and chosen maximum current for DIII-D.Here we present new 3D electromagnetic response calculations performed with the ThinCurr model [18] which include the response of the coils and the surrounding conducting structure to a prescribed current rampdown.The structures are modeled as a set of 2D elements coupled by a series of circuit equations.The 3D models for  3 shows that the spectra do not differ much at very early times when the non-linearly growing modes dominate, but differ at late when the applied fields dominate.The ideal wall boundary condition is maintained for modes n = 7 − 10 even in the resistive wall case.
the conducting structures are able to capture major features like portholes, and varying resistivity for the elements allow different material thicknesses to be represented.The SPARC and DIII-D models include over 15 000 and 38 000 elements, respectively.The plasma is represented as a set of filaments in order to model the response to a linear current rampdown.In the SPARC modeling with ThinCurr, resistance is added to the coil in order to limit the maximum coil current to 350 kA (well below the previously predicted 590 kA) due to engineering constraints (specifically a limit on tolerable j×B forces).In the DIII-D modeling, a maximum current just below 100 kA is predicted, slightly less than the minimum chosen for previous calculations.In both cases, a nonlinear fit to the coil current vs. plasma current replaces the previous linear assumption.The results of the REMC coil current response for a ThinCurr simulation of each device is shown in figure 6.The time response and 3D fields from these calculations are used for NIMROD modeling in the next section.But, as we will describe below, we do not use the calculated fields at the FW directly; rather, we use the fields at the LCFS as a basis to determine the appropriate boundary condition.

Simulations with extrapolated normal fields at the VV
Whether ideal or resistive, the conducting wall boundary condition is more appropriately applied at the location of the VV and not the FW.The fact that the coil is interior to the VV prevents the coils fields from being used as boundary conditions at the VV directly.As an alternate strategy, we calculate a set of normal magnetic field boundary conditions at the VV that would give approximately the same field distribution at the LCFS as that produced by the coils, if there were instead a vacuum in between.In all calculations, only the normal magnetic fields at the boundary are prescribed, with the fields interior determined by the governing equations.The LCFS is chosen as the boundary at which to match the applied fields as the natural boundary between the vacuum and plasma regions.In these initial calculations, the method for extrapolating the fields to the VV is to choose a set of virtual toroidal circular coils outside of the VV.For each Fourier component individually, we solve for the currents (necessarily having non-zero divergence for n > 0) in the virtual coils that would produce the desired (ThinCurr) fields on the LCFS, then calculate the normal fields produced by these currents at the VV boundary.A result of this calculation for SPARC is shown in figure 7.Although the virtual coil approach is able to produce field distributions at the LCFS that closely match the desired fields, a more direct method of calculating the extrapolated fields using the current potential [29] at the outer boundary may be employed instead in future calculations.

SPARC simulations with VV boundary
For SPARC we compare two simulations, both beginning with the fields calculated using ThinCurr.In one case, the ThinCurr fields are applied directly at the original FW-shaped boundary, and in the other case, the extrapolated fields at the VV boundary (from the ThinCurr fields at the LCFS) are used.Both simulations model only the CQ phase of the disruption, neglecting the TQ-induced MHD.In figure 8 we compare the evolution of the magnetic energy spectra for the two cases as well as the decay of the plasma current.Note that by necessity, the vacuum fields in the region between the FW and VV will be larger in amplitude in the VV boundary case, so that that globally integrated magnetic energy is expected to be larger in the VV case even for similar amplitudes in the plasma region.However, qualitative differences can be seen in the time evolution of the energy spectra.Whereas the nonlinear growth and saturation of the n = 1 − 10 modes in the FW boundary case has one relatively broad peak at 0.5 ms, the VV-boundary simulation has an earlier sharp peak in the mode amplitudes, followed by a second, broader peak.Although a current spike is not expected when the TQ phase is not included, the sharp peak in the mode amplitudes is accompanied by a significant inflection in the plasma current trace in the VV-boundary simulation, associated redistribution of the current profile.
The RE test particle population used in SPARC resembles that in DIII-D, except that the range of initial energies is from 1 to 50 Mev, because at higher magnetic field, the upper limit of RE energies that can be confined before curvature drift becomes comparable to device minor radius increases.The initial sharp spike in the n = 1 − 10 magnetic energy corresponds to a rapid and total loss of RE test-particles, whereas a more gradual loss of the entire population occurs with the ThinCurr fields applied at the FW boundary (see figure 9).The early rapid loss in the VV case corresponds to a loss of good flux surfaces throughout the volume, but this is followed by a brief reformation of good flux surfaces just after the initial sharp peak in the magnetic energy.The second broader peak corresponds to a second loss of confinement that persists for a longer time.Although this 'multiple crash' behavior was seen in prior simulations that included the TQ (where the first crash corresponded to the TQ-induced MHD and not to the REMC perturbations directly), the reformation and subsequent re-stochastization of the flux surfaces has not been seen in CQ-only simulations with the small FW-shaped boundary.A secondary test-population was not launched, but it would be expected that any new population relaunched before 1 ms would again be lost, corresponding to more than 50% of the CQ.

DIII-D simulations with VV boundary
The grid geometry used for the VV-boundary simulations in DIII-D (figure 10) does not extend the simulation domain as far outside the FW as in the SPARC case.However, the simulations that were performed with the IWL plasma shape by definition had the FW-shaped boundary extremely close the LCFS, so we may expect significant changes in the behavior with this constraint relaxed.Two simulations with the VV grid are performed, beginning with the IWL equilibrium and considering only the CQ phase of the disruption.In one case the vacuum fields that were used for the original IWL simulations (with 100 kA maximum current) in [19] are extrapolated from the LCFS to the VV boundary, for direct connection with those results.In the other case, the coil fields calculated by ThinCurr at the LCFS are extrapolated to the VV.
The overall RE test-particle loss curves for the two simulations (figure 10) are initially different but by 1.4 ms end up overlapping very closely.In both cases, the RE losses are continuous and exceed 99%, which is quite dissimilar from the previous IWL simulations in which the loss fraction plateaued at around 90% in every case.We can also see the loss is step-wise, which is more evident in the lower curves, showing the derivatives of the loss curves and an oscillatory loss rate (which on average exceeds the avalanche growth rate).The greater separation of rational surfaces in the IWL configuration, which was responsible for the insensitivity to coil current in the original results, also tends to produce this series of discreet loss events as successive rational surfaces are destroyed.The earlier losses seen with the ThinCurr fields appears to be an effect of a relatively higher n = 2 amplitude compared to the n = 1 when the effect of all the conducting structures from ThinCurr is included.As the time when the edge is initially stochasticized in the simulations with ThinCurr fields, n = 2 islands can be observed, where only wellseparated n = 1 islands are present with the vacuum fields (figure 11).
For only the case with the ThinCurr fields, the RE loss curve is also plotted broken down by energy range.The testparticles are seeded with uniformly distributed energies from 1-50 MeV and the individual energies do not evolve.It can be seen that at early times in the simulation the loss rate for mid-range energy REs (10-30 Mev) is greater than REs at higher or lower energy.This produces a significant valley in the distribution function close to 20 MeV (figures 10(d) and (e)), although this is not representative of an expected distribution function in an experiment, since in reality all REs would have to accelerate through the 20 MeV bottleneck.Preferential loss of mid-energy REs is the combined result of two physical effects.First, the higher energy REs have orbits that are shifted outward in major radius due to curvature drift ).The relativistic factor γ causes the extent of this shift to increase with increasing RE energy.The outward shift increases the probability of the a particle being lost to the outer wall.Second, the diffusive loss on stochastic field decreases with increasing RE energy [30].This occurs because the shifted orbits tend to average over smaller magnetic perturbations.In this simulation, only REs in the lowest energy bin are predominantly lost to the inner wall.In the higher two energy bins, the REs are lost to the outer wall due to a combination of the shift in major radius and diffusion in minor radius.The two effects, with opposite energy dependence, add in such a way that the region in which REs with 20 MeV can remain confined has a smaller radius than REs of either higher or lower energy.

Conclusions
REMC simulations for SPARC and DIII-D have been updated to include more realistic boundary conditions in three respects: first, the ideal wall model was replaced with a resistive wall; second, the ThinCurr code was used to compute more realistic boundary fields based on the response of the coil and surrounding conducting structures to a prescribed current quench; third, the location of the boundary was moved outward to the approximate location of the VV, beyond the internal REMC.In order to enable this, an equivalent set of 3D fields had to be calculated to apply at the VV boundary that would produce approximately the same field distribution at the LCFS assuming a vacuum between the two.The first of these changes was explored independently, with ideal and resistive wall simulations compared directly in the FW shape geometry.A final set of simulations incorporated these two changes in addition to the resistive wall, with each change isolated for one of the two machines.In SPARC, a comparison of the ThinCurr fields for the FW boundary and extrapolated to the VV was performed, whereas for DIII-D, the comparison was between two cases with fields extrapolated to the VV: the ThinCurr fields and the vacuum fields.
The comparison of ideal and resistive wall simulations for DIII-D showed an effect of the boundary condition on the reformation of flux surfaces in the core after a TQ in a LSN diverted geometry.In a certain respect, it may be surprising to see this kind of effect on the core stochasticity from the distant boundary, but the effect is not primarily mediated by the stability of internal modes.Rather, the resistive wall allows the applied spectrum of external modes to evolve through redistribution of the fields toroidally and poloidally.This in turn changes the degree of resonance of the applied spectrum with the plasma.In the particular case of the DIII-D simulation, the reduced 2/2 and enhanced 3/2 increased resonance near the core.While flattening of the mode spectrum can not be assumed to increase resonance in every case, it would be expected to allow some effect of the coil over a broader range of q-profiles.SPARC also showed enhanced deconfinement of REs with a resistive wall, but in that case the effect was due to delayed reformation of flux surfaces in the core.The delay was correlated with a more gradual increase in q on the axis after the TQ, allowing resonance with the coil to be maintained longer.
While the use of a boundary condition with the appropriate wall resistance will always be more accurate than the ideal wall approximation, the extent to which this is true will depend on the relationship between the disruption timescales and the wall time.For both of these machines, the wall time is short enough to be relevant on the disruption timescale.Locating the conducting wall the appropriate distance from the plasma should be important regardless of wall time.
For the DIII-D IWL equilibrium, extending the boundary to the VV made a significant difference in the loss fraction.With the FW shaped boundary, losses never exceeded ∼90%, while both the vacuum fields and the ThinCurr fields extrapolated to the VV produced continuous losses in excess of 99%.A non-monotonic energy dependence was also observed in the RE losses, where REs with energy near 20 MeV were lost faster than those at higher or lower energy.This occurs due to the combination of outward shifted orbits (increasing the propensity to strike the outer wall relative to lower energy REs) and diffusive transport on stochastic fields (reduced for even higher energy REs), resulting in a smaller confined region for REs in this energy range than higher or lower energies.
In the case of SPARC, simulations with the ThinCurr fields for both the FW and VV geometries showed a complete loss of the RE test-population, but the test particles were lost earlier and more rapidly with the VV boundary.Further, the VV boundary simulation showed a 'multiple-crash' behavior, in which the flux surfaces were destroyed by an initial MHD event and began to reform before a second MHD event fully stochastized the fields once again.In other simulations that included the TQ, we have seen an initial MHD event due to the TQ and a second due to the coil, but this case included only the CQ, and the multiple-crash behavior was not seen in previous CQ-only simulations.
With only RE test-particles, a notable lack in the present model is the effect of any RE currents on the MHD evolution.The goal is to prevent the formation of any appreciable contribution of the REs to the total current, but only so long as that is achieved does the model remain accurate in this respect.If a significant RE current begins to form in any region of the plasma, this would effect the overall current profile and therefore the resonant interaction with the coil.Since we have seen that rising q on axis decreases the effectiveness of the coil, a growing RE current localized to the core could have a positive feedback effect of bringing the central q back down below 2, increasing resonance with the coil, and allowing the nascent RE population to be deconfined each time it begins to form (sawtooth-like behavior).However, such behavior is speculative and a model including RE currents is needed to determine whether the level of RE current required to produce this effect would already be dangerously large.
These simulations demonstrate that boundary conditions play a critical role in accurately modeling REMC performance, with each of the model improvements here resulting in an improvement in the expected performance.The goal of these improvements is to have a model that can be carefully validated once REMC data on SPARC and DIII-D becomes available, and then used as a predictive tool for future tokamak designs.

Figure 1 .
Figure 1.The geometry of the REMCs designed for SPARC and DIII-D, along with the vacuum vessel model used in the ThinCurr code (discussed in part 4).Reproduced from [18].© 2023 The Author(s).Published by IOP Publishing Ltd on behalf of the IAEA.All rights reserved CC BY 4.0.

Figure 2 .
Figure 2.With a resistive wall (lower row), the SPARC disruption has a more gentle increase in q 0 after the Ip-spike, compared to the ideal wall simulation (upper row).In both simulations, some flux surfaces reform in the core shortly after the TQ, but as the REMC grows during the CQ, the core flux surfaces are again destroyed when q 0 remains below 2.

Figure 3 .
Figure 3.In DIII-D simulations, RE confinement is significantly reduced with a resistive wall boundary condition, even more so with the REMC than for a disruption with no coil.(a) Thermal energy and plasma current for two disruption simulations with no REMC (dashed) and two with an REMC (solid).The REMC has an effect on the global evolution, whereas the ideal (blue) or resistive (red) wall has very little effect.(b) Magnetic energy in the n = 1 (thicker lines) and n = 2 (thinner lines) for the four simulations.In the TQ, n = 2 initially dominates over n = 1.(c) Confined REs vs time for the four simulations, where the initial population is launched at t = 0, then lost during the TQ.A new population is launch 0.5 ms after the end of the TQ in each case, which is earlier with the REMC included.

Figure 4 .
Figure 4.A comparison of the late time (3.0 ms) toroidal magnetic energy spectrum for the ideal and resistive wall DIII-D REMC simulations, showing transfer of the applied n = 1 energy to higher n components when the normal fields at the wall can redistribute.The comparison of just n = 1 and n = 2 in figure3shows that the spectra do not differ much at very early times when the non-linearly growing modes dominate, but differ at late when the applied fields dominate.The ideal wall boundary condition is maintained for modes n = 7 − 10 even in the resistive wall case.

Figure 5 .
Figure 5.A comparison of flux surface reformation in DIII-D afterthe TQ shows that the appearance of a good confinement region near the axis is delayed with a resistive wall.The field lines in the Poicaré plots are color coded by the number of toroidal transits before striking the wall, where 20 000 transits is the maximum integration length.

Figure 6 .
Figure 6.REMC current response for ThinCurr simulations of SPARC (left) and DIII-D (right).In each simulation, a filament representation of the plasma current is self-similarly ramped-down linearly (green lines).The SPARC coil has added resistance to keep the REMC current (blue lines) limited to 350 kA.In each case, a nonlinear curve is fit to the REMC current as a function of the plasma current (dashed red lines) during the rampdown phase.This relationship between REMC amplitude and plasma current is used in the NIMROD modeling.

Figure 7 .
Figure 7. (Left side) Normal magnetic fields are applied at a smoothed version (mangenta) of the SPARC vacuum vessel (VV) shape (grey dashed) to achieve approximately the same normal field distribution at the LCFS (blue) as what is calculated directly by ThinCurr.Black dashed line is the SPARC FW shape.(Right) Real (solid) and imaginary (dashed) parts of the n = 1 normal magnetic field distribution for each point along the discretized boundaries, beginning at outer midplane and proceeding clockwise.The fields calculated by ThinCurr at the LCFS (blue) are used to solve for equivalent fields to apply at the VV (magenta).Given these applied fields, the actual fields resulting at the LCFS are shown in red, which are meant to closely match the ThinCurr fields.

Figure 8 .
Figure 8. (Upper) Toroidal mode spectrum of δB/B ( √ Wn/W 0 ) for SPARC simulation with fields from ThinCurr applied at the FW boundary.(Middle) δB/B spectrum from SPARC simulation with ThinCurr fields applied at the VV boundary.(Lower) Decay of the total current for the two simulations, where black is FW boundary and magenta is the VV boundary.

Figure 9 .
Figure 9. (Upper) Confined test population of REs vs time for the two SPARC simulations with ThinCurr fields, where black is for the FW shaped boundary and magenta is with the fields extrapolated to the VV boundary.(Lower) Magnetic field line Poincaré plots are shown for several times (as indicated) for the VV boundary simulation.Field lines are colored by the number of toroidal transits before striking the boundary.After the loss of all REs a brief reformation of flux surfaces in the core occurs, followed by a secondary loss of confinement.Since new REs can be produced at any time during a real CQ, confinement after the initial loss of the test population remains important.

Figure 10 .
Figure 10.(a) The shape of the LCFS for the DIII-D limited equilibrium (blue), the shape of the DIII-D FW (black-dashed) and the shape of the larger VV grid boundary outside the FW (magenta).(b) RE test-particle loss results with the field extrapolated to the VV grid boundary.The heavy cyan line is the simulation with the vacuum fields (used in prior simulations) extrapolated to the new boundary, while the heavy blue line is simulation using the extrapolated ThinCurr fields.The black, red, and yellow curves are the results from the ThinCurr case broken down by energy, where black is less than 10 MeV, red is 10-30 MeV and yellow is more than 30 MeV.The population is uniformly distributed in energy from 1-50 MeV.(c) Total RE loss rates (derivatives of the upper curves) compared to the avalanche growth rate (dashed).(d) and (e) Histogram of number of confined REs in energy bins from 1-50 MeV at time 1.0 and 1.6 ms respectively.A preferential loss of midrange energy REs is initially seen.

Figure 11 .
Figure 11.Magnetic field line Poincaré plots in straight-field-line coordinates for two IWL DIII-D simulations, at the time when the edge begins to stochasticize for the simulation with field calculated by ThinCurr (left).The simulation using the vacuum fields (right), which are dominated by the n = 1 component, has only isolated n = 1 islands at the same time.The ThinCurr fields have an n = 2 component almost as large as the n = 1.