Runaway electron dynamics in ITER disruptions with shattered pellet injections

This study systematically explores the parameter space of disruption mitigation through shattered pellet injection in ITER with a focus on runaway electron (RE) dynamics, using the disruption modeling tool Dream. The physics fidelity is considerably increased compared to previous studies, by e.g. using realistic magnetic geometry, resistive wall configuration, thermal quench onset criteria, as well as including additional effects, such as ion transport and enhanced RE transport during the thermal quench. The work aims to provide a fairly comprehensive coverage of experimentally feasible scenarios, considering plasmas representative of both non-activated and high-performance DT operation, different thermal quench onset criteria and transport levels, a wide range of hydrogen and neon quantities injected in one or two stages, and pellets with various characteristic shard sizes. Using a staggered injection scheme, with a pure hydrogen injection preceding a mixed hydrogen-neon injection, we find injection parameters leading to acceptable RE currents in all investigated discharges without activated runaway sources. Dividing the injection into two stages is found to significantly enhance the assimilation and minimize RE generation due to the hot-tail mechanism. However, while a staggered injection outperforms a single stage injection also in cases with radioactive RE sources, no cases with acceptable RE currents are found for a DT-plasma with a 15MA plasma current.


Introduction
Disruptions -off-normal events when the plasma energy is rapidly lost -and associated generation of highly energetic runaway electron (RE) beams represent an outstanding challenge of the tokamak concept for magnetic confinement.
An effective disruption mitigation system in a tokamak reactor should limit the exposure of the wall to localised heat losses during the thermal quench (TQ), and to the impact of RE beams, as well as avoid excessive forces on structural elements.In an ITER plasma at 15 MA plasma current, avoiding excessive electromagnetic forces requires the current quench (CQ) time to be kept in the range 50-150 ms, and the RE current should be kept below 150 kA [1].
The currently envisaged mitigation method is to inject a massive amount of material, primarily consisting of a mixture of neon (Ne) and hydrogen (H), when an emerging disruption is detected, in order to better control the plasma cooling and energy dissipation.
In ITER, the baseline strategy for delivering the material is through shattered pellet injection (SPI) [1].However, finding operation parameters in reactor-scale devices that simultaneously avoid RE formation and mitigate heat and electromagnetic loads remains an open question.
Several previous studies of integrated simulations of the full disruption event have been made with a focus on the CQ evolution, based on simplified geometry and/or assumptions of a prescribed material deposition and temperature drop [2,3,4].Although earlier results [2] indicated that the RE generation could be suppressed by a sufficiently large H injection, later results suggested that the recombination occurring at very high H quantities counteracts the damping of the runaway generation, thereby boosting the avalanche RE generation [3,4], leading to several MAs of RE currents in all studied cases of ITER-like DT plasmas.
An initial study of simulations of SPI-mitigated disruptions in ITER was presented in [5], using a selfconsistent model for the injection and temperature drop for a given evolution of the magnetic perturbation amplitude.The main objective of this paper was to investigate the effect of separating the injection into two phases, as previously considered in [6].The resulting two-stage cooling was found to effectively thermalize the tail of the electron distribution function before the onset of the RE generation, strongly reducing the RE seed due to the hot-tail mechanism.This indicated a possibility to obtain low RE currents in plasmas without tritium decay and Compton scattering of photons from the radioactive wall (referred to as activated sources).Nonetheless, large RE currents were still found for DT plasmas.However, these simulations studied a rather limited set of scenarios and still included several simplifications, such as circular flux surface cross sections, a fixed TQ onset time (independent of the evolution of the plasma parameters), and neglected ion transport once the injected material was ablated.
In this study we expand the parameter space as well as the range and complexity of the scenarios compared to the previous studies, while also improving the physics fidelity.The simulations are performed with the disruption simulation tool Dream [7].The Dream code is designed to include much of the relevant physics during a disruption while keeping the complexity sufficiently low (e.g., by resolving one spatial dimension) to make it feasible to explore a large parameter space, or even allowing for optimization [4].It includes an SPI ion source based on an analytical ablation rate [5], time dependent rate equations for ionisation and radiation, as well as a self-consistent evolution of the induced electric field.
Dream covers a wide range of the computational fidelity-affordability trade-off concerning RE physics, up to a linearized bounce-averaged relativistic Fokker-Planck treatment.Here, we use a computationally inexpensive fluid RE model, allowing the exploration of a large parameter space.Further details about the model and the simulation settings are covered in section 2 and the simulation results are described in section 3. The results and limitations of the model are discussed in section 4, before the paper is concluded in section 5.

Simulation Settings
The simulations are based on four reference scenarios produced by the CORSICA [8] workflow, which are commonly used within the ITER Organisation.These scenarios are referred to as DTHmode24, a highperformance H-mode D-T plasma (pure, and of equal isotope concentrations) with a 15 MA plasma current; H26, a hydrogen L-mode plasma with a 15 MA plasma current; He56, a helium H-mode plasma with a 7.5 MA plasma current; and H123, a hydrogen H-mode plasma with a 5 MA plasma current.The flux surface geometry and initial profiles for the studied cases are illustrated in Fig. 1.Moreover, we assume a resistive wall time of t wall = 0.5 s and a wall radius of r wall = 2.833 m, chosen to match the poloidal field energy within the vacuum vessel (which is the toroidally closed conductor closest to the plasma) in JOREK [9] simulations, corrected for contributions from the poloidal field coils.

SPI parameters
The hydrogen isotope of pellets used in our simulations is deuterium (D) which is expected to produce very similar results to H, the isotope planned to be employed in ITER.The reason for this choice is that the pellet ablation model we use is developed for D. The pellet injection speed is assumed to be v p = 500 m/s, and the fragment velocity dispersion is uniform within v p ±∆v, with ∆v/v p = 0.4.The injection is modelled in a poloidal plane within a spreading angle of α = 10 • , with the shards diverging from the point {R, Z} = {8.568,0.6855} m (with the origin at the horizontal midplane at the center of the torus), as shown in Fig. 1d.We take into account a B-field correction of the ablation, representing the diamagnetic shielding of the pellet [10].
The SPI fragment sizes follow the distribution proposed by Parks [11], and we consider three different characteristic shard sizes, corresponding to 487 (default shard size), 5185 (small shard size) and 68 (large shard size) shards for a standard ITER pellet (consisting of 1.85•10 24 atoms), as also done in previous INDEX simulations [12].As the Ne concentration is varied, the D quantities are adjusted to keep the total number of atoms in the pellet constant.We consider both single and multiple pellet injections with up to four pellets, injected in either one or two stages.
We consider two Ne concentrations for each scenario, chosen to roughly span the desired ohmic CQ time range of ∼ 50 − 150 ms, to the extent possible (see related discussion about Fig. 3a in section 3.2.1).The ohmic CQ time is defined here as t CQ = [t(I Ohm = 0.2I p0 ) − t(I Ohm = 0.8I p0 )] /0.6, where I Ohm is the total ohmic current, and I p0 is the initial total plasma current.For each case t CQ is evaluated based on an additional simulation with runaway generation disabled but otherwise identical settings.The choice of Ne concentrations is based on initial tests on a baseline case, using a single stage injection, the default shard size and a TQ of 3 ms duration with the "late" onset criterion, (see section 2.2).The chosen Ne quantities are N Ne,inj = 2 • 10 23 and N Ne,inj = 1.83•10 24 for the H26 scenario; N Ne,inj = 2.5 • 10 22 and N Ne,inj = 1.5 • 10 24 for the DTHmode24 scenario; N Ne,inj = 5 • 10 21 and N Ne,inj = 1.83 • 10 24 for the He56 scenario; and N Ne,inj = 10 21 and N Ne,inj = 1.83•10 24 for the H123 scenario.Additionally, for every plasma scenario considered, we perform simulations using an injected Ne content of N Ne,inj = 5 • 10 22 , to assist comparison between scenarios, as well as to connect the results to previous work with INDEX [12].
The cloud of pure hydrogenic pellets is prone to build an excess pressure and thus drift significant distances towards the low field side [13,14,15].Although a computationally efficient model to selfconsistently capture this effect has been developed [16], it is not yet available in Dream.Instead, to emulate this effect in two stage injections modelled here, the material deposition of the first set of pellets (with no neon content) is shifted one radial grid cell (with a width of about 20 cm) outward of the shard position ‡.Then the pellet shard plume experiences a considerably reduced dilution cooling generated by the material deposition, leading to an outward-shifted, sometimes quite edge-localised, deposition profile.To bracket the range of drift effects, for comparison we also perform similar simulations with local deposition.
The ablated material is initially deposited into the neutral charge state, and the ionization is then evolved using time dependent rate equations.The ionization and recombination rates, as well as the corresponding radiation rates, are taken from the ADAS database [17] for Ne and the AMJUEL database § for D. The latter accounts for opacity to Lyman radiation, which has previously been found to be important at the high injected D densities studied here [5].

Thermal quench
We consider two thermal quench onset criteria.According to the early onset criterion, the TQ is triggered when a Ne containing shard with a horizontal velocity of v p reaches the q = 2 surface.In the late onset criterion the TQ is triggered when any point inside of the q = 2 flux surface reaches a temperature of 10 eV; in case of pellets with low neon concentration this criterion is usually, although not always, satisfied significantly later than the early one -hence the name.
These criteria are motivated by the point that the resistive decay of the current profile is sufficiently fast at 10 eV to likely trigger a magnetohydrodynamic (MHD) instability around the q = 2 flux surface on a millisecond time scale.The temperature drop to 10 eV, however, might occur much earlier locally (i.e., somewhere within the flux surface) than on flux surface average.In particular, such a local cooling might happen in the vicinity of Ne containing shards, which, in an extreme case, might be sufficient to trigger an MHD instability as the shards pass the q = 2 flux surface.This scenario motivates the early onset criterion, while the late onset criterion corresponds to the opposite extreme case, where the temperature variation within the flux surface remains negligible.These two criteria are thus intended to bracket the range of possible TQ onset delays.
Once the TQ is triggered, spatially and temporally homogeneous electron heat diffusion and runaway electron particle diffusion are activated, and remain active for a period of t TQ = 1 ms or 3 ms in different scenarios.These time scales are in the range indicated by extrapolations of empirical data from present day devices, and can be considered as typical values used to assess thermal loads in an ITER TQ [18].The diffusion coefficients are parameterized by a relative magnetic perturbation amplitude, δB/B, using a Rechester-Rosenbluth-type model [7,19] for the heat transport, and the model by Svensson et al. [20] for RE transport.The values of δB/B are chosen such that the temperature at the magnetic axis would have time to drop below 200 eV in the time t TQ if transport would be the only loss mechanism (note though, that it typically reaches much lower values due to radiative energy losses).At lower temperatures, Ohmic heating can balance the heat losses due to transport, preventing the temperature from dropping much further due to transport alone.The resulting values of δB/B are tabulated in Appendix A, and are of the order of 0.1-1% -comparable to those extracted from an MHD simulation of an SPI triggered ITER TQ, in [21].
Ion particle transport (of all ion species and charge states) is activated at the same time as the other transport channels, with initial diffusion and advection coefficients of D ion = 4000 m 2 /s and A ion = −2000 m/s, which then decay exponentially over the characteristic time of τ = 0.5 ms.A similar form of the time evolution of the ion transport coefficients was used in [22] to reproduce experimental data in the ASDEX Upgrade tokamak with the ASTRA code.The chosen values of D ion and A ion give diffusion and advection time scales (r 2 minor /D ion and r minor /A ion , respectively, where r minor is the minor radius of the plasma) in the ms range.Together they result in a transport of a substantial amount of material to the plasma core over ∼ 0.1 ms, as expected from 3D MHD simulations [21], and the balance between the advection and the diffusion gives a moderately peaked final profile (r minor /L n ≈ −r minor A ion /D ion ≈ 1 where L n is the characteristic density length scale).The decay time constant τ = 0.5 ms ensures that the ion transport is significant on the same time scale as the other transport processes are active.

Runaway electrons
The Dream simulations use the fully fluid runaway model.In simulations of the H26, He56 and H123 scenarios, only non-activated runaway seed generation mechanisms are included: the Dreicer generation uses a neural network trained on kinetic simulations∥ [23], and the hot-tail model used is described in Appendix C in [7].The avalanche multiplication rate is described in [24], and it accounts for both partial screening and magnetic trapping effects.The Svensson model [20] is used to turn the momentum p and pitch ξ = p ∥ /p dependent diffusion coefficients into diffusion coefficients of the runaway fluid.Here the assumed p-dependence is p/(1 + p 2 ), that reproduces the ∝ v low energy behavior, and emulates a reduction of the runaway transport due to finite Larmor radius and orbit width effects at high p.The ξ-dependence corresponds to the ∝ |v ∥ |-dependence of the Rechester-Rosenbluth type diffusion coefficients.In addition, the DTHmode24 case also includes RE generation by activated sources.
According to [2], the photon flux from the wall is expected to drop by a factor of ∼ 1/1000 immediately after the neutron bombardment of the wall is stopped (i.e., when the fusion reaction ceases).To account for this in the calculation of the Compton seed, we start the simulations with the nominal photon flux on ITER given in [2], and then reduce it at the end of the TQ (taken to be the end of the transport event) by a factor of 1/1000.The ITER nominal photon flux taken from [2] was calculated assuming a beryllium wall design.

Simulation results
The difficulty of avoiding a large RE formation is strongly dependent on the initial plasma current, due to the maximum RE avalanche gain being exponentially sensitive to it.We therefore study the scenarios with a 15 MA plasma current separately in section 3.1 and 3.2, before turning to the cases with reduced plasma current in section 3.3.

Representative example
We start by considering one of the baseline cases (late TQ onset criterion and a 3 ms TQ) with the H26 scenario, which serves as a representative example of a single stage injection.Fig. 2 shows the time evolution of the most important quantities for the case when the neon quantity is adjusted (10.8%) to achieve a CQ time of 100 ms.The total hydrogenic ion density increases along the trajectory of the fastest shards (lower solid green line), as shown in Fig. 2a.After the ions are rapidly mixed across the entire radius in the beginning of the transport event (indicated by the lower dotted line), the density remains steady for the rest of the simulation.Note that in single pellet injection cases the density evolution of neon is very similar, within the multiplicative factor of the injected Ne/D ratio.We also note that the assimilated fraction of the injected material in this case is f assim = 5.18%, which corresponds to 1.04•10 22 assimilated neon atoms and 8.55 • 10 22 hydrogen atoms.
In Fig. 2b, we can see the temperature drop below the keV-range after the first shards reach the plasma, and once the temperature locally goes below 10 eV just inside the q = 2 flux surface, the transport event is initiated.This leads to a more homogeneous further cooling of the plasma, ending at T e ∼ 10 eV.The initially created RE seed that follows the first cooling gets mostly lost in the transport event, then the RE current recovers and increases slowly and reaches macroscopic values around t = 40 ms, as shown with the dash-dotted line of Fig. 2c.Note that the RE current at the end of the transport event is tiny, ∼ 10 −10 A (outside the scale of figure 2c), and it mostly stems from a hot-tail seed.A small, but non-negligible, Dreicer seed generation remains active during the earlier parts of the CQ, generating a representative seed of ∼ 10 −8 A, but the seed still remains very small.
The representative RE current is still significant, 5.31 MA, due to the very strong avalanche.We introduced two quantities; the representative seed is the total RE current at the end of the TQ (that is not lost due to transport) plus the seed generated after the TQ; and the representative RE current is RE current when 95% of the total current is carried by REs.Henceforth these are understood when we refer to a scalar RE current or seed without further qualifications.
An important general feature -which is more pronounced in the well performing cases -is that the temperature drop is usually divided into two separate stages, despite all material being injected in a single pellet.This can effectively suppress the hot-tail seed generation, as this enables the hot-tail to slow down and thermalize with the bulk at an intermediate temperature before a strong electric field is induced after the second, radiatively dominated, cooling phase.
In addition, in a late TQ the significantly higher density increase, which is allowed by a higher temperature, also assists thermalization.
The separation of the cooling stages has also been observed in simulations with the INDEX and JOREK codes [25], and is explained by the temperature dependence of the various heat loss mechanisms.Transport as a heat loss channel becomes inefficient in the 100 eV range and below, while radiative cooling due to the presence of neon has a local minimum around 200 eV and then rises sharply towards lower temperatures (see e.g.Fig. 6 in [5]).During the injection phase, the temperature is also strongly affected by dilution.Due to the 1D nature of the employed model, the flux surface averaged dilution cooling occurs instantaneously after deposition, thus the effect of dilution on the temperature is determined directly by the ablation rate.The strong scaling of the ablation rate with temperature makes this process dominant in the early stage of the injection (before the onset of the transport event), while this process also becomes inefficient in the 100 eV range and below.
In case of a late thermal quench, the cooling starts with a dilution down to a few hundred eV, and then it slows down due to the low radiated power.Eventually, the temperature drops below 100 eV in some parts of the plasma where cooling accelerates again due to radiation.
As the temperature drops below 10 eV at some radial location, the transport event is triggered.Transport-induced cooling, although remaining moderate in this temperature range, helps drive the rest of the plasma below 100 eV, initiating a global radiative collapse.This global radiative collapse is further assisted by the particle transport, increasing dilution, as well as the radiating neon content, in regions with low or no direct deposition of injected material.
The separation of the cooling phases is typically less pronounced in cases with an early TQ, and also in cases with large neon injections, typically resulting in larger RE currents.In those cases, the earlier cooling also reduces the ablation, and hence the frictional drag on energetic electrons, which further enhances the RE generation.Indeed, the scenario shown in Fig. 2a-c when simulated with an early TQ gave a representative RE current of 7.72 MA.When aiming for a CQ time of 50 ms (which in this case means a pellet consisting of 99% neon, see section 2.1), we obtained a representative RE current of 11.7 MA.
On the other hand, both the separation of the cooling phases and the collisional damping of RE generation are more pronounced with a staggered injection.In those cases, the lack of radiating neon during the first injection stage allows for a higher temperature in the vicinity of the shards, resulting in a higher ablation, and also delays the onset of the transport event and the radiative collapse with both the early and late TQ onset criterion.An example of such an evolution is shown in Fig. 2d-f for a similar case to that in Fig. 2a-c but with the neon-containing pellet preceded by a full deuterium pellet by 5 ms.Notably, in this case the RE seed becomes negligible, and hence no RE current is generated.
The drift of the ablated material can lead to a shift in the deposition profile from the core towards the edge.This shift can be much larger than the drift length itself, since as the dilution cooling caused by the shards is displaced behind the plume of shards, the higher temperature in the vicinity of the shards greatly enhances their ablation.Thus the shards may ablate away rapidly, before reaching the core.
In the example shown in Fig. 2d, where the shift between the ablation location and the deposition location is only ∼ 20 cm, the outward shift of the deposition profile is moderate, which is due to the moderate temperature of the background plasma of the H26 scenario.In contrast, in the ∼ 4 times hotter DTHmode24 plasma we sometimes observe deposition profiles peaked as far out as r ∼ 1.5 m, for the same deposition shift of ∼ 20 cm (not shown here).The deposition may in some cases be so localized, that the corresponding dilution cooling may be sufficient to meet our late TQ trigger condition, even before the neon doped shards arrive into the plasma.This is the case in the example shown in Fig. 2e, where the transport onset (lower dotted line) appears at an earlier time than the second injection (upper pair of solid lines).
In spite of these, potentially dramatic, differences in the early evolution of the disruption between local and shifted material deposition, the effects of the shifted deposition are reduced by the transport event.Namely, the density profile is strongly flattened by the fast ion mixing during the transport event.As a result, the CQ evolutions with and without the shift included are rather similar.We should keep in mind though, that our results rely on an effective ion mixing; the efficiency of this process is a source of uncertainty.Moreover, if the pellet cloud drift would be comparable to the plasma minor radius it could cause a significant amount of material to be completely ejected from the plasma, which might have a significant effect on the CQ evolution.It is worth noting that within our model, assuming that no such material ejection takes place, the shift also leads to a somewhat increased amount of assimilated D in the plasma, as will be discussed in connection to Fig. 4.

Correlation among the explored parameters with a 15 MA plasma current
In this section we study correlations between the various settings and figures of merit covered in this study, by plotting the various settings and metrics against each other, displaying pairs of quantities with observable trends of interest.A list of the simulation settings included is provided in Appendix B. Section 3.2.1 focuses on quantities related to the plasma current evolution, while section 3.2.2focuses on quantities related to the assimilation of the injected material.
3.2.1.Current evolution Fig. 3 shows two correlation images, with different markers for the H26 and the DTHmode24 simulations, with single stage and staggered injection and, for the DTHmode24 scenario, with and without activated sources.In Fig. 3a we can see that there is a slight trend of reduced RE current with increasing current quench time.Such correlation is expected as a longer t CQ tends to be accompanied by lower induced electric fields (during both the TQ and the CQ) to which all the runaway generation processes, in particular the Dreicer and hot-tail mechanisms, are sensitive.
Notably, choosing the Ne content aiming to span the t CQ = 50 − 150 ms range in the baseline cases keeps t CQ within or in the vicinity of the desired range also when the parameters are varied.However, we see that the H26 (blue and green circles) and the activated DTHmode24 (black and red squares) cases occupy different regions in the I RE − t CQ -space.The H26 cases appear at higher t CQ values and span a wider range in I RE , extending from ∼ 12 MA down to negligible RE currents.The DTHmode24 cases where the activated sources are artificially removed (orange) also have negligible RE currents, while all DTHmode24 cases with activated sources included (black and red) have RE currents above ∼ 4 MA.
It was not possible to span the entire 50 − 150 ms CQ time range for each scenario and injection scheme.In the H26 scenario, even a full Ne pellet was insufficient to reach the lower end, and for the DTHmode24 scenario, reaching the upper end would typically require too low Ne content, for which the TQ is incomplete (except for a couple of cases with a staggered injection and shifted deposition).In these incomplete TQ cases a significant part of the plasma is ohmically re-heated to hundreds of eV following the TQ, which results in an intolerably long CQ.
The difference between the H26 and DTHmode24 scenarios depends mostly on the difference in initial thermal energy.In H26, the relatively low thermal energy content limits the amount of material ablated before the temperature becomes too low to allow any significant further ablation (see Fig. 4 and the corresponding discussion in section 3.2.2).
It is thus not possible to assimilate a sufficient amount of material to achieve a post-TQ temperature corresponding to the lower part of the desired CQ range.In contrast, in DTHmode24, the relatively high thermal energy content increases the amount of assimilated Ne, resulting in a low post-TQ temperature, while increasing the cooling required to initiate a radiative collapse.In this scenario the CQ times thus populate the lower part of the desired range, while decreasing the injected quantities in an attempt to reach the upper part often leads to an incomplete

TQ.
Sometimes the ohmic current profiles resulting from an incomplete TQ are more localized, and are expected to be prone to instabilities.While worth further investigation in the future, the stability of the current channels is not monitored in our simulations.Instead, we consider all incomplete TQ cases as undesirable and we choose the Ne concentrations to avoid them.Incomplete TQs nevertheless occurred in a few cases even if the TQ was complete in the baseline cases examined when choosing the Ne concentrations; these are the cases plotted with open symbols at far right side of the scale.Note, however, that although these cases have a very long CQ time according to our definition (i.e, measured in identical simulations without RE current), some of these actually develop significant RE currents which abort the CQ preemptively, thus their CQ time is not representative.
Negligible RE currents occur mainly because the RE seeds are extremely small in those cases.In fact, in most such cases the RE currents calculated in the simulation correspond to much less than a single electron (marked by the vertical line in Fig. 3b) ¶, indicating that with high probability no REs would be generated in those cases at all.
Importantly, even a tiny seed can be amplified to ¶ Note, that Dream, being a continuum code, is agnostic to the discreteness of charges, thus it can produce currents much smaller than that corresponding to a single relativistic electron.However, the cases studied here are not ambiguous in that they either produce much smaller or much larger seed currents -corresponding either to no REs, or a sufficiently large RE population that the continuum treatment is justified.Note, that in Fig. 3b the points shown left of the "single RE in ITER" line are in fact located far to the left of this limit.
several MAs of final RE current.This is illustrated in Fig. 3b, showing an essentially logarithmic dependence of the final RE current on the representative seed, enabling even a seed as small as 10 −8 A to be amplified to ∼ 5 MA.To lead the eye, the trends found in the DTHmode24 and H26 cases are shown with dashed lines, and while the values and the slopes are comparable between the two plasmas, we find that the same seed typically yields a higher final RE current in the H26 plasma.This must be caused by a higher avalanche rate, which is consistent with the typically poorer hydrogen assimilation in H26 caused by the lower initial temperature (see Fig. 4), which leads to lower effective critical electric fields.In addition, while the H26 cases extend to very small seeds (and even to zero seed REs, outside the scale of Fig. 3b), in the DTHmode24 cases the seed cannot be reduced below the minimum representative Compton RE seed of ∼ 10 −3 A. As a result, we find no RE currents below ∼ 4 MA in DTHmode24.

Assimilation
The lower CQ times seen in figure 3a for the DTHmode24 cases can be linked to higher ablation figures caused by the significantly higher initial thermal energy, as mentioned in section 3.2.1.This is reflected in the assimilated Ne and D quantities shown in Fig. 4a, where the DTHmode24 cases typically occupy higher values in both Ne and D than H26 (with an exception for a couple of the staggered DTHmode24 cases with low Ne assimilation).
We also see that there is a negative, approximately exponential, correlation between the assimilated Ne and D, with different slopes for different groups of points.The negative trend occurs because a smaller assimilated Ne quantity reduces the radiative cooling,  so that the temperature remains higher in the vicinity of the shards, allowing for a faster ablation and thus a higher assimilated D quantity.Conversely, a larger amount of D gives a stronger dilution cooling, which slows down the ablation and thus reduces the amount of assimilated Ne.
For the single stage injections (blue for H26, black and orange for DTHmode24), the different slopes correspond to the cases with early and late TQ.For the staggered cases (which were only performed with a late TQ), the division occurs between the cases with and without the shift of the deposition included.Another clear trend is that the staggered cases typically have larger amounts of assimilated D and, consequently, lower amounts of assimilated Ne.
The RE current is rather strongly negatively correlated with the amount of assimilated D, as shown in Fig. 4b.This is partly due to the negative correlation between the amount of assimilated Ne and D, as smaller Ne quantities lead to later TQ onsets, longer separation between the cooling phases and lower electric fields, all of which reduce the RE generation.A higher D quantity also has the direct effect of increasing the collisional damping.As a result, both the staggered H26 cases and the DTHmode24 cases without activated sources all have negligible RE currents.There are also two single stage H26 cases with negligible RE currents, one of which has the lowest amount of injected Ne among the H26 cases (N Ne,inj = 5 • 10 22 ), and one of which uses three pellets.Both of them used TQ conditions favourable for RE avoidance, i.e. a late TQ and t TQ = 3 ms, leading to a relatively large ablation, long separation between the cooling phases and small hot-tail seed.
Increasing the number of shards (i.e. using the small shard size) can also increase the D assimilation, but the effect is rather moderate (∼ 10%).Moreover, it also increases the assimilated Ne, which can even lead to a moderate increase in the RE current.However, the staggered cases appear to more robustly suppress the RE formation.
At the highest assimilated D quantities obtained with the DTHmode24 scenario, the decreasing trend of I RE with N D,assim flattens out, and a hint of an emerging increasing trend can be seen.The reason for this is that large enough D quantities can make the hydrogen species recombine, while there is still a significant ohmic current left in the plasma that can be converted to REs.At high electric fields, this decreases the collisional damping while retaining the number of available targets for the RE avalanche, resulting in an increase in the RE current, as discussed in [3].

Cases with reduced plasma current
We now turn our attention to the scenarios with lower initial plasma current (He56 and H123).As the maximum avalanche gain scales exponentially with the plasma current, we now expect significantly smaller conversion rates for a given seed compared to the 15 MA scenarios studied above.This is indeed confirmed by the simulation results, although it still seems important to use a staggered injection to robustly avoid a significant RE formation for all realistic TQ conditions.
Figure 5 shows representative current evolutions following a single stage and a staggered injection with unfavourable (early and short) TQ conditions, for a) the He56 scenario and b) the H123 scenario.For the He56 scenario, we inject N Ne = 5 • 10 21 Ne atoms for both the single stage and staggered injection.In the staggered injection case with the H123 scenario, we inject N Ne = 1.83 • 10 24 Ne atoms.For the single stage injection, we decrease the Ne quantity and increase the D quantity to assess whether it is feasible to avoid a significant RE formation even with a single stage injection if one would be allowed to exceed the upper CQ time limit.In the case displayed we use 3 pellets with a total of N Ne = 10 21 Ne atoms.
We see in Fig. 5a) that a disruption with unfavourable TQ conditions may result in several MAs of RE current even with a plasma current of 7.5 MA, if it is only mitigated by a single stage injection; in this case the representative RE current is 2.48 MA.Similarly, in figure 5b) we see that the RE current could not be reduced much below 1 MA with a single stage injection even with a 5 MA initial plasma current, despite exceeding the upper CQ time limit; the displayed case has a representative RE current of 0.77 MA.
The representative RE seed is close to 3 A in both single stage injection cases, demonstrating the potential severity of the RE avalanche also at these reduced plasma currents in ITER.The mitigation is also somewhat impeded by a rather low assimilation (∼ 10%) for the single stage injections.This is similar to the corresponding cases with the H26 scenario, despite the ∼ 2 times higher temperature, as the density is significantly lower, making the plasma easier to cool.
In the staggered cases, however, we see no significant RE formation, regardless of the Ne quantity in the second stage (as long as the CQ remains complete).The same applies to single stage injections with favourable TQ conditions, except for extreme cases with very large Ne quantities.We also note that the RE conversion is considerably smaller for a given seed than in the 15 MA scenarios; at a representative seed of 3 A, we find a conversion of about 50% for the 15 MA scenarios (see figure 3b), while the single stage injection cases in figure 5 have conversions of 33% for the He56 scenario and 15% for the H123 scenario.

Cases with injection arriving after the TQ
Finally we consider cases where the injection is triggered too late, so that it arrives after the TQ transport event; specifically after a 50 ms delay.All the late injection cases with standard injection settings lead to re-heating and incomplete current quench, as illustrated in Fig. 6, which is based on the H26 scenario.Here we use conditions that would favor avoiding a re-heating, namely inject a 99% neon pellet and employ a remnant electron heat transport after the TQ corresponding to δB/B = 4•10 −4 and 2•10 −3 (note that the latter is comparable with δB/B during a TQ).We see that after the TQ, the temperature increases due to ohmic heating to T e ≲ 1 keV with a remnant δB/B = 4 • 10 −4 , and T e ≲ 200 eV with a remnant δB/B = 2 • 10 −3 .When the injection arrives, the temperature drops to the 10−100 eV range, but a large part of the plasma re-heats again to a temperature of a few hundred eV.In other cases we used a fixed remnant  heat diffusivity of 1 m 2 /s, and obtained qualitatively similar results.The reason for the re-heating is the extremely low assimilation rate of the pellet, owing to the low plasma temperature (< 1 keV).Thus, increasing the remnant heat diffusivity does not help, as it just decreases the temperature further, and along with it the assimilation rate (it is 0.4% for δB/B = 4 • 10 −4 and 0.2% for δB/B = 2 • 10 −3 ).Similarly, moving the time of the injection closer to the transport event, when the plasma temperature has had less time to regrow, further aggravates the re-heating problem.In conclusion, it is very important for the utility of SPI that the pellets do not arrive late with respect to the TQ.
The re-heating can however, in some cases, be counteracted by reducing the shard size, and thereby increasing the assimilation.This was possible in an He56 scenario with t TQ = 3 ms, a 99% Ne pellet injection and a fixed remnant heat diffusivity of 1 m 2 /s, where we increased the number of shards from 487 to 5185 (corresponding to the standard and small shard sizes as named in [12]).This increased the assimilation rate from 0.19% to 0.31% that, in this case, was sufficient to eliminate the re-heating issue.Nevertheless, reducing the shard size to the range considered here does not appear to be a robust method to avoid re-heating, although the assimilation might be further enhanced by employing a more extreme shattering.To ameliorate this problem ITER plans to employ three upper injectors dedicated to post-TQ injection, aiming at injecting mainly gas by using a very large shattering angle.Assessing gas penetration is however out of scope of the present study.
It is also important to mention that the case where re-heating was avoided produced an almost full runaway conversion; as no injected material had reached the plasma at the time of the TQ, the collisional drag remained rather low, making the plasma very prone to RE seed generation.Here we have however assumed the TQ duration to remain in the 1-3 ms considered for mitigated disruptions, while the TQ duration in ITER in the absence of impurities might be substantially longer, up to several tens of milliseconds [26].This could significantly reduce the RE seed generation, and also expand the time range during which the plasma is hot enough to significantly ablate the pellet, allowing for a later pellet arrival.

Discussion
The results presented here are based on modeling choices which are, in many ways, conservative, and it is thus possible that, in particular, the RE currents calculated here are overestimated.One reason for this is that Dream does not account for the vertical motion of the plasma and the resulting shrinking of the closed field line region.As the REs will very quickly get lost to the wall on the open field lines, the vertical motion acts to decrease the volume with unhindered runaway formation.
Although a detailed study of the effect of a vertical displacement event is outside the scope of the paper, one can get an indication of the impact of this effect by varying other aspects of the boundary condition at the conducting vessel wall, with the aim to emulate the evolution of the poloidal flux available for RE generation on a given flux surface.In particular, setting a perfectly conducting wall very close to the plasma, thus forcing the poloidal flux at the edge to remain constant, could reduce the poloidal flux variation and the resulting avalanche gain considerably.Such a boundary condition would also make the simulated poloidal flux variation on a given flux surface similar to the actual variation during the time the flux surface is closed, if the plasma moves in such a way that the poloidal flux at the last closed flux surface remains approximately constant.This is expected to be the case if the resistive time scale of the plasma wall is long compared to the CQ time, which will be true for mitigated disruptions in ITER, and is also observed in recent ITER simulations with the JOREK code [27].
To illustrate the range of possible outcomes in case the poloidal flux available for runaway generation is modified by some process, figure 7 shows a comparison between the current evolutions obtained with different wall boundary conditions.The baseline case with an effective conducting wall radius of r wall = 2.833 m and resistive wall time of t wall = 500 ms is compared to a case with r wall = r minor (where r minor is the plasma minor radius) and t wall = 500 ms, and one with Time evolution of the plasma current and RE current using different settings for the wall boundary condition, for the staggered DTHmode24 case using 1 pure D pellet for the first injection and a pellet with 1.35% Ne for the second injection (target CQ time 100 ms for the corresponding baseline case).r wall = r minor and t wall = ∞.The comparison is made for the D-T case with the lowest RE current included in figure 4b (excluding the cases without the shifted deposition for pure D injections).This case uses a staggered injection with 1 pure D pellet for the first injection and a Ne concentration of 1.35% in the second injection.
Reducing the wall radius to r wall = r minor is found to reduce the RE current by several MAs.When also setting t wall = ∞, the RE current becomes as low as ∼ 1 MA.While using t wall = 500 ms increases the RE current to ∼ 2 MA compared to the t wall = ∞ case.The RE current also keeps increasing during the plateau phase as additional poloidal flux can diffuse through the wall.Conversely, with t wall = ∞, a clear RE decay is observed.It may be however, that the behavior observed here during the plateau phase is not realized in an experiment, since the flux surfaces may completely disintegrate before the RE plateau is fully formed.
It should be noted that the reduction might not be quite as large in reality; for instance, the magnetic energy carried by the escaping REs will not immediately be lost from the system, but will remain in the form of a halo current, and the electric field associated with this halo current might diffuse into the confined plasma region where it can contribute to the RE generation.Nevertheless, this exercise indicates that the effect of the plasma motion, and in general modified wall boundary conditions, may result in a major reduction in the generated RE current, motivating more detailed future studies on this front.
An additional conservative modelling choice is that no RE transport losses due to magnetic perturbations are included during the CQ.While the overall perturbation level is expected to be much smaller than during the TQ, a δB/B as low as ∼ 10 −4 may be sufficient to have a notable impact on the RE generation and decay [20].An even larger uncertainty stems from the unknown time of reformation of the outermost flux surfaces after the TQ.The RE current may also be reduced when accounted for by a kinetic calculation, as the fluid RE generation rates, in particular for hot-tail generation, have been shown to overestimate the RE generation [28].
There are also modelling choices which may lead to an underestimation of the RE generation, although perhaps less significant compared to the overestimation effects mentioned above.For instance, in this study we have not accounted for the rapid current density profile relaxation and the corresponding current spike typically taking place during the TQ.This effect has been shown to increase the RE generation in cases prone to RE generation in the outer part of the plasma [29].However, a shift of the RE current density profile towards the edge would also make the REs more prone to be lost to the wall during a vertical displacement event; the result of the combination of these competing effects remains to be investigated.
Another uncertainty regarding the RE generation is that the γ-photon flux and spectrum was calculated for a first wall made of beryllium, while ITER currently plans to have a first wall made of tungsten.This might alter the γ-photon flux, and consequently the Compton scattering RE seed generation.However, as observed in figure 3b), the final RE current (if macroscopic) is only logarithmically sensitive to the RE seed; approximately 0.4-0.9MA change in the final RE current per order of magnitude change in the seed.Thus the γ-photon flux must change by orders of magnitude before it would lead to a major difference in the final RE current.Moreover, a recent assessment shows that the difference in the prompt γ-photon flux and spectrum between a tungsten and a beryllium first wall is not significant [30].While the corresponding differences in the decay γ-photon flux and spectrum present during the CQ remain to be assessed, current estimates indicate that the change of the first wall material is not likely to have a major impact on the RE generation.
Finally, there are several sources of uncertainty which are not directly related to the RE generation.One of these concerns the details of the transport event.We have attempted to cover the plausible range of onset criteria and magnetic perturbation levels, but recent calculations indicate that the TQ time scale can be longer than assumed here.This could reduce the hot-tail seed considerably, possibly allowing for a complete RE avoidance in a wider range of scenarios with a single stage injection without activated RE sources.Moreover, the ion transport coefficients are a significant source of uncertainty.While our values are successfully chosen to match the mixing time scale found in 3D MHD simulations, this argument is not sufficient to accurately constrain the involved parameters, and thus both the time scale of the ion mixing and the resulting radial profile might vary considerably.This could have a particular importance for the staggered injection cases with a shifted deposition, which are dependent on the ion mixing for material deposition in the core.In addition, it is worth emphasising that the shift of the deposition has here been assumed to be moderate, while there might occur shifts comparable to the plasma minor radius [16], which could eject substantial amounts of material from the plasma and thereby reduce the assimilation.

Conclusion
We have explored a wide range of feasible disruption mitigation SPI settings and conditions for four ITER scenarios with a focus on RE avoidance: one L-mode hydrogen plasma with a 15 MA plasma current (H26); one H-mode helium plasma with a 7.5 MA plasma current (He56); one H-mode hydrogen plasma with a 5 MA plasma current (H123), all without activation of the wall; and one DT H-mode scenario with a 15 MA plasma current, including activated runaways sources (DTHmode24).The disruption mitigation performance was evaluated based on the CQ time, which should be in a given range -50 − 150 ms in 15 MA discharges -and the RE current, which should be lower than 150 kA.We observe a wide range of runaway conversions from 80% at the largest Ne and lowest D quantities and unfavourable TQ conditions, down to negligible values for low Ne quantities, large D quantities and favourable TQ conditions.
The most important conclusions that can be drawn from this study are the following: (i) The two-stage injection is an effective means to suppress the runaway current.It suppresses the hot-tail seed generation by leaving time for the tail of the distribution to thermalize, and it reduces the avalanche rate by improving hydrogen assimilation.
In particular, with two-stage injection the runaway current can be completely eliminated in non-activated discharges even at 15 MA initial plasma current.
(ii) While two-stage injection also yields the best performing cases in activated scenarios, there the irreducible nuclear seeds -active even after the transport losses of the TQ -yield a floor of ≈ 4.5 MA in the final RE current (corresponding to a mA-level seed).
(iii) Even in the reduced initial current scenarios it is possible to get multi-MA runaway currents, especially for unfavorable (early and short) TQ conditions.However, in these cases the twostage injection can robustly eliminate the runaway current.(iv) Both concerning t CQ and the final RE current, it is important that the pellet injections arrive before the TQ.Otherwise, the poor ablation in the much cooler background plasma leads to insufficient material assimilation and a corresponding ohmic re-heating of the plasma, along with an intolerably long CQ time.Although the poor assimilation may be ameliorated by post-TQ gas injection and extreme pellet shattering, this will likely not be sufficient to compensate for the lack of collisional damping of the RE seed generation during the TQ.Thus, unexpectedly high runaway currents may be generated.(v) The current quench time can be kept within the desired range in all plasma scenarios considered, with appropriate injected quantities.(vi) The final RE current, if it is macroscopically large, is logarithmically sensitive to the effective runaway seed current (i.e., the RE current surviving the TQ losses, plus the seed generated afterwards).Megaampere-scale currents may arise from a single runaway electron in the avalanchedominated regime of ITER.(vii) The external electric field boundary conditions have a major effect on the final runaway current.This includes the location and the resistive timescale of the wall, as well as effective changes in these boundary conditions that arise due to a scraping off of the current channel during vertical displacement events.(viii) A negative correlation with an exponential sensitivity of the assimilated neon quantity to the assimilated hydrogen quantity is observed.
Staggering tends to yield higher hydrogen and smaller assimilated neon quantities.The overall assimilation figures are significantly higher in the hot DTHmode24 scenario than in the H26 scenario.(ix) A high assimilated hydrogen content is beneficial in most cases; it is a requisite to obtain scenarios with negligible runaway current.However at the highest assimilated values (few times 10 24 assimilated atoms), which are achievable only in DTHmode24, we observe the previously predicted increase of runaway current caused by hydrogen recombination.
The results presented here are based on conservative modelling choices, and they should prompt con-tinued experimental validation of Dream, as well as further studies accounting for additional effects.For instance the available poloidal magnetic flux can potentially be reduced, and runaway losses enhanced, during a vertical displacement event.In addition, the possibility of RE instabilities and remnant magnetic perturbations after the thermal quench may also reduce the RE current.

Figure 2 .
Figure 2. Plasma parameter evolutions for a baseline H26 case with a target CQ time of 100 ms.a-c) single stage injection, and d-f) staggered injection with a shifted deposition for the first injection.a) and d) Hydrogen density (injected+background).b)and e) Electron temperature.c) and f) Time evolution of the total plasma current (dashed), and its ohmic (solid) and runaway (dashdotted) components.In a-b) and d-e) solid green lines indicate the trajectory of the fastest and slowest shards for each injection, the dashed line marks the q = 2 flux surface, and the temporal boundaries of the transport event are shown with dotted lines.
DT single st., no activated seed H single stage, long CQ H staggered, local dep.DT staggered, local dep.H staggered, shifted dep.DT staggered, shifted dep.

Figure 3 .
Figure 3. Correlations between figures of merit related to the plasma current evolution, varying all the simulation parameters within the ranges described in section 2. a) Representative RE current vs. CQ time.Cases with incomplete CQ are shown with open symbols at t CQ = 160 ms.b) Representative RE current (linear) vs. representative RE seed (logarithmic).Blue circles -single stage H26, light green circles -staggered H26 w. shifted deposition, dark green circles -staggered H26 w. local deposition, black squaressingle stage DTHmode24, light red squares -staggered DTHmode24 w. shifted deposition, dark red squares -staggered DTHmode24 w. local deposition, orange diamonds -single stage DTHmode24 w.o.activated RE seed.In panel b) the seed corresponding to a single relativistic electron is marked by a vertical line; all points to the left of that line are far to the left of the scale plotted range; these points are plotted there to illustrate that the final RE current in those cases is negligible.The typical range of activated seed generated after the TQ is indicated by the gray shaded region.
DT single st., no activated seed H single stage, long CQ H staggered, local dep.DT staggered, local dep.H staggered, shifted dep.DT staggered, shifted dep.

Figure 4 .
Figure 4. Correlations between figures of merit involving assimilation of the injected material.a) Number of assimilated Ne (logarithmic) vs. D atoms (linear).b) Representative RE current (linear) vs. number of assimilated D atoms (logarithmic).Blue circles -single stage H26, light green circles -staggered H26 w. shifted deposition, dark green circles -staggered H26 w. local deposition, black squares -single stage DTHmode24, light red squares -staggered DTHmode24 w. shifted deposition, dark red squares -staggered DTHmode24 w. local deposition, orange diamonds -single stage DTHmode24 w.o.activated RE seed.

Figure 5 .
Figure 5. Current evolutions following a disruption with unfavourable TQ conditions mitigated with a single stage injection and a staggered injection, for a) the He56 scenario and b) the H123 scenario.For the He56 scenario, the injected Ne quantity is N Ne = 5 • 10 21 atoms with both the single stage and the staggered injection.For the H123 scenario, the injected Ne quantity is N Ne = 1.83 • 10 24 atoms for the staggered injection, and N Ne = 10 21 atoms for the single stage injection, in the latter case injected with 3 pellets.

r
Figure 7.Time evolution of the plasma current and RE current using different settings for the wall boundary condition, for the staggered DTHmode24 case using 1 pure D pellet for the first injection and a pellet with 1.35% Ne for the second injection (target CQ time 100 ms for the corresponding baseline case).