Effect of Two-Stage Shattered Pellet Injection on Tokamak Disruptions

An effective disruption mitigation system in a tokamak reactor should limit the exposure of the wall to localized heat losses and to the impact of high current runaway electron beams, and avoid excessive forces on the structure. We evaluate with respect to these aspects a two-stage deuterium-neon shattered pellet injection in an ITER-like plasma, using simulations with the DREAM framework [M. Hoppe et al (2021) Comp. Phys. Commun. 268, 108098]. To minimize the obtained runaway currents an optimal range of injected deuterium quantities is found. This range is sensitive to the opacity of the plasma to Lyman radiation, which affects the ionization degree of deuterium, and thus avalanche runaway generation. The two-stage injection scheme, where dilution cooling is produced by deuterium before a radiative thermal quench caused by neon, reduces both the hot-tail seed and the localized transported heat load on the wall. However, during nuclear operation, additional runaway seed sources from the activated wall and tritium make it difficult to reach tolerably low runaway currents.


Introduction
The sudden loss of the energy confined in fusion plasmas during off-normal events, called disruptions, presents one of the most severe threats to the future of fusion energy based on the tokamak design. An efficient disruption mitigation system is therefore of utmost importance for future high-current devices such as ITER. The potentially greatest threat to be mitigated is posed by large currents carried by highly energetic runaway electrons, whose number increases exponentially during the runaway avalanche [1] from any small seed population, and which may cause severe damage upon wall impact. The disruption mitigation system must also ensure sufficiently homogeneous deposition of the thermal energy on the plasma-facing components, and avoid excessive forces on the machine resulting from currents flowing in the surrounding structures.
The mitigation method currently envisaged is the injection of massive quantities of material when an emerging disruption is detected, aiming to better control the plasma cooling and energy dissipation. Conventionally, the material is delivered as a gas puff from a pressurized vault [2]. This technique, while comparatively simple, has a number of disadvantages. The injected gas ionizes rapidly when exposed to the initially hot plasma, and so becomes tied to the magnetic field, substantially slowing the spread of the gas through the plasma.
Moreover, the change in plasma profiles during the gas injection can accelerate the growth of plasma instabilities, allowing the disruption to progress before injected material has reached all parts of the plasma.
Another approach, that can provide better penetration, is to inject material in the form of a solid, cryogenic pellet. The exposure to the plasma causes the pellet to ablate and deposit its content along its trajectory. The ablation can be made more efficient by shattering the pellet into smaller shards before it enters the plasma. This shattered pellet injection (SPI) technique has been chosen as the baseline for the disruption mitigation system at ITER [3].
The details of the design and the operation parameters for the disruption mitigation system in reactor-scale devices remain an open question. A complexity arises as different aspects of the mitigation generate potentially conflicting requirements. For example, radiative dissipation of the thermal energy favours large injected quantities, which may result in an unacceptably fast quench time of the plasma current, due to the high plasma resistivity at low temperature. In addition, rapid cooling increases the hot-tail runaway seed [4,5].
A method suggested recently to circumvent these issues is to divide the injection into two stages [6].
First, injection of a large amount of pure deuterium is used to cool the plasma through dilution, by a factor 10-100, without significantly perturbing the magnetic field configuration or radiating any substantial amount of thermal energy. A few milliseconds later a smaller neon injection follows, which radiatively dissipates the thermal energy. The division of the cooling into two steps gives the hot electrons (in the tail of the electron distribution) time to equilibrate at an intermediate temperature before the runaway generation is initiated, potentially suppressing the hottail runaway generation. In addition, if the plasma perturbation destroying the plasma confinement is initiated at a lower temperature, the slower thermal motion reduces the thermal energy transport. Instead, a larger fraction of the thermal energy can be dissipated through radiation, reducing the danger of localised wall hot spots.
The results presented in reference [6] indicate that it is indeed possible to substantially cool the plasma through dilution by a deuterium shattered pellet injection, without destroying the plasma confinement. However the outcome and optimisation of this twostage injection scheme, particularly the runaway generation, have not yet been studied in detail.
Recent simulations of plasma shutdown scenarios, focused on the dynamics during the current quench, indicate that material injection can lead to high runaway currents in ITER [7,8]. While a massive increase in the electron density appeared as a promising route to suppress the runaway avalanche [7], the substantial recombination in such scenarios results in a net enhancement of the final runaway current [8]. These studies used a simplified instantaneous impurity deposition with a flat radial profile and prescribed thermal quench dynamics.
The hot-tail runaway seed generated during the thermal quench varies over many orders of magnitude, depending on the injected densities and seed losses arising from perturbations of the magnetic field [9,10]. Regions of parameter space with sufficiently small hottail seeds have been found previously, which are not expected to result in an excessive final runaway current despite the strong avalanche gain [10]. However, these studies did not model the consistent build-up of the injected density following SPI and the subsequent current quench dynamics.
In this work we investigate the effect of SPI on runaway generation in a tokamak disruption. We use the recently developed numerical tool Dream (Disruption Runaway Electron Analysis Model) [11], which is capable of self-consistently calculating the time evolution of the background plasma properties, the electron momentum distribution and the runaway current during a mitigated tokamak disruption. In particular, we simulate two-stage SPI in an ITERlike setting, and assess how the disruption mitigation performance depends on the amount of injected deuterium and neon. This assessment is based on a quantification of the radiated power, current quench time and the maximum runaway current, which must remain within defined limits for machine protection.
We find that the maximum runaway current is a non-monotonic function of the injected deuterium density, with an optimal range identified, that shows a strong sensitivity to the opacity of the plasma to Lyman radiation. The two-stage injection is found to be effective in reducing both the hot-tail seed and the transported energy losses. However, during nuclear operation, additional runaway seed sources, tritium decay and Compton scattering, make it difficult to reach tolerably low runaway currents.
The rest of this paper is structured as follows: in section 2, we briefly introduce Dream and detail the SPI model that is used in this work. Section 3 describes the effect of two-stage SPI with different injection parameters, with an emphasis on how the disruption dynamics is influenced by dividing the injection in two stages, the difference between nuclear and non-nuclear operation, and the importance of opacity to Lyman radiation. The influence of the model assumptions on the results presented is discussed further in section 4, and conclusions are summarised in section 5.

Shattered pellet injection and plasma dynamics
The evolution of the plasma parameters and runaway electrons during the disruption mitigation are modelled here with the 1D-2P numerical tool Dream [11]. Dream was extended in this work to include the capability to model SPI based on a Neutral Gas Shielding (NGS) model [12].
The employed modelling accounts for the plasma response to the SPI, transport of thermal energy by magnetic perturbations and the kinetic evolution of the electron distribution function, capturing the hottail as well as the Dreicer runaway generation. In the non-nuclear phase of tokamak operation, we find that the runaway seed is dominated by hot-tail generation. In the nuclear phase of operation, additional runaway seeds are provided by tritium decay and Compton scattering of gamma photons from the radioactive wall. These are modelled as quasi-stationary sources feeding particles directly into the runaway population [13]. After the thermal quench, when the hot-tail mechanism is no longer active, the Dreicer mechanism is also modelled in a similar fluid-like fashion [11].
The amplification of the runaway seed by the avalanche mechanism is modelled accounting for effects of partial screening [14]. This introduces a significant dependence of the runaway avalanche on collisionalradiative processes and the resulting distribution of ionization states. As we will show, the response to a large deuterium SPI is significantly influenced by the opacity of the plasma to Lyman radiation. For convenience, an overview of the particle and energy balance, and self-consistent electric field evolution, implemented in Dream are given in Appendix A.
An SPI starts with a pellet being accelerated and then shattered, resulting in a plume of pellet shards of different sizes and velocities entering the plasma. As the shards travel through the plasma, they act as a set of localised sources of particles. The corresponding density increment-here calculated assuming instantaneous homogenisation over the flux surfaces-depends on the ablation rate of the shards, the area of the flux surfaces, and a kernel function defining the radial spread of the deposited material around the shards.
The technical aspects needed to describe the SPI scenario are presented below in more detail. These include the size and velocity distributions of the shattered pellet shards, the explicit form of the ablation rate, and how that is translated to a density source.

Shattering
We assume the pellets are shattered into N s approximately spherical shards, with sizes described by an equivalent radius r p,k , with k = 1, ..., N s . The shard radii are drawn randomly from the distribution with probability density P (r p,k ) = k 2 p r p,k K 0 (k p r p,k ), where K 0 is the zeroth modified Bessel function of the second kind, k p = N inj /6π 2 n p N s −1/3 , n p is the number density of the solid pellet material, and N inj is the total number of injected atoms [15]. This distribution of shard sizes has recently been used in several other SPI studies [16,17,18]. Once shattered, the shards are assumed to travel with constant velocities v p,k in a poloidal plane, starting at the shattering point (x 0 , y 0 ), as illustrated in figure 1; with the origin of the (x, y)-coordinate system taken at the plasma center. The speeds v p,k , and angles α p,k with respect to the horizontal plane, are chosen from uniform random distributions within v p ± ∆v p and ±∆α p , respectively. The parameters x 0 , y 0 , v p , ∆v p and ±∆α p , as well as N inj and N s , are considered controllable free parameters. In reality, the control of these parameters is achieved by adjusting the impact speed and angle of the pellet on the shattering surface, likely resulting in correlations between them.

Ablation
We characterise the ablation rate by the time derivatives of the shard radii,ṙ p,k . The expression used forṙ p,k is based on a version [12] of the NGS model [19] that allows the pellet material to have both hydrogenic and noble gas components. Expressed in terms of the unidirectional incident heat flux q in carried by the bulk plasma electrons and their effective energy E in , this model giveṡ Here, the normalising radius, heat flux and effective energy are r p0 = 2 mm, q 0 = n 0 2T 3 0 /(πm e ) and E 0 = 2T 0 , respectively, with the representative temperature and density T 0 = 2000 eV and n 0 = 10 20 m −3 . The solid mass density of the pellet is denoted ρ dens . The dependence of the ablation rate on the deuterium-neon composition is accounted for by the factor λ(X) = [27.0837 + tan (1.48709X)]/1000 kg/s, where X = N D2 /(N D2 + N Ne ) is the deuterium fraction, N D2 is the number of deuterium molecules (thus the number of deuterium atoms is N D = 2N D2 ) and N Ne is the number of neon atoms in the pellet.
The heat flux and effective energy are calculated from a general electron momentum distribution function, f , according to and E in = 2 n free m e c 2 (γ − 1)f dp.
We note that equation (2) was derived assuming a Maxwellian electron momentum distribution, with temperature T M . However, it may be assumed to be sufficiently accurate for the small deviations from a Maxwellian present in the early stages of the disruption while the shards are still ablating, that is, before a substantial runaway acceleration has occurred ‡. The total free electron density is n free = f dp, c is the speed of light, and γ is the Lorentz factor. The factor 1/4 in equation (3) converts the isotropic heat flux to the average unidirectional heat flux facing the pellet shards, and is strictly valid for a Maxwellian distribution [21]. The normalisation of E in is chosen such that E in reduces to 2T M for completely Maxwellian electrons, which is equal to the ratio of the unidirectional heat flux and the unidirectional particle flux.

Material deposition
The ablated material quickly ionizes, becoming confined by the magnetic field to the flux surfaces near the pellet shard position where the ablation took place. The homogenization and equilibration of the ablated material is approximated here to take place instantaneously, an assumption also made in other recent SPI studies [6,17,18,22]. The impact of this assumption is discussed later in this paper.
The homogenized ion density increase on the flux surface with radius r is given by where the pellet molar mass is M, and N A is the Avogadro number. A fraction f ij of the ablated pellet material appears in the charge state i of ion species j, with a corresponding density n ij . The radial distribution of the deposited material, in terms of density increase, is described by the factor H(r, ρ p,k ) = h(r, ρ p,k )/A fls (r), where A fls = 4π 2 rR 0 is the area of the flux surface at radius r and h(r, ρ p,k )dr describes the fraction of the material deposited at a radius between r and r + dr ablated from a pellet located at minor radius ρ p,k . The width of this deposition kernel is physically determined by transport processes radially dispersing the ablated material. Some previous studies have used a Gaussian deposition kernel h ∝ exp [(r − ρ p,k ) 2 /r 2 cld ] [6,16], with a "cloud width" of r cld that has a lower bound of ∼ 1 cm [23], the width of the flux tube channeling the pellet material in the vicinity of the shard. After verifying that r cld has only a weak quantitative impact on the final density profile, we opted for a delta function deposition kernel, h = δ(r − ρ p,k ) for numerical convenience, translating to a uniform distribution over the distance travelled by the shard during one time step. ‡ In the presence of non-negligible runaway populations different physics mechanisms dominate [20].
To avoid the need to resolve the extremely rapid ionization dynamics of the ablated material, we assume the deposited material appears in its equilibrium distribution of charge states, at the local density and temperature. With φ j denoting the particle fraction of the pellet material consisting of species j, the quantity f ij appearing in equation (5) becomes f ij = φ j n eq ij /n tot,j , where the equilibrium distribution of charge states is calculated according to i n eq ij = n tot,j , i = 0, 1, ...Z.
The total density of species j, with atomic number Z, is denoted by n tot,j , and I ij (T M , n M ) and R ij (T M , n M ) are the ionization and recombination rates, respectively, obtained from the OpenADAS database [24]. The energy required for the initial ionization is accounted for by a loss term (A.3) in the energy transport equation (A.1).

Plasma response to SPI
Dream calculates the time evolution of the pellet and plasma parameters self-consistently, within a model summarised as follows. The pellet shards follow straight lines as described by equation (1), and the evolution of the shard sizes is governed by the ablation rate, equation (2). The resulting density increase is given by equation (5), and the resulting cooling is calculated by the heat transport equation (A.1) with the additional loss term given in equation (A.3).
Once the pellet material is deposited, the evolution of n ij is governed by ionization and recombination according to the time dependent rate equations (A.5). The electric field evolution resulting from the rapid change in the conductivity is given by the induction/diffusion process (A.7). The momentum distribution of the electrons is resolved using the "fully kinetic" mode described in [11].

Simulation of two-stage shattered pellet injection
In this section we evaluate the disruption mitigation performance of the two-stage deuterium-neon injection scheme in an ITER-like plasma. With this scheme, the plasma is first cooled by dilution down to the 100 − 1000 eV range by a deuterium injection, without significantly affecting the thermal energy density (and hence the pressure).
A few ms later, neon is injected, with the aim of radiating away most of the thermal energy, and producing a current quench with an acceptable timescale. The temporal separation between the injections allows the tail of the electron distribution to thermalize, thereby reducing hot-tail runaway generation. Another advantage is that the ratio of the transported and radiative energy losses decreases if the thermal quench is triggered from a lower temperature, reducing the risk of damage in plasma facing components due to excessive localized heat loads. These aspects are therefore emphasised in our investigation.

Injection and plasma parameters
We search for suitable injection parameters for ITERlike plasmas with initial temperature profile T M (r) = T c (1 − 0.99(r/a) 2 ), central temperature T c = 20 keV, and a radially constant initial density n M = 10 20 m −3 . The profiles used here have previously been used to study massive material injections in ITERlike scenarios [7,8], assuming flat deposition profiles. For the initial ion composition, we consider both a pure deuterium plasma and an even mix of deuterium and tritium, the latter corresponding to the nuclear operation phase. The plasma minor radius is a = 2 m, the radius of the perfectly conducting wall is b = 2.15 m and the major radius of the tokamak is R 0 = 6.2 m. The initial plasma current is I p = 15 MA, with radial profile given by A fast reaction time of the disruption mitigation system favours high pellet speeds. We therefore consider the fastest injection speeds expected. These are v p,D = 800 m/s for deuterium pellets and v p,Ne = 200 m/s for neon pellets [3]. For the distribution of pellet speeds we assume ∆v p = 0.2 v p and we set the divergence angle ∆α p = 20 • , for both deuterium and neon pellets. This angle primarily affects how many shards pass through the innermost flux surfaces, and an increased divergence can shift the deposited density profile from the core further out in the plasma. While this may be of interest when fine-tuning the density profile, we are able to achieve satisfactory density profiles with this fixed value. The current quench dynamics is expected to be rather insensitive to the details of the density profile for a given number of injected particles, as long as core penetration is achieved [8]. We also fix the position at which the pellets are shattered to the tokamak wall on the horizontal mid-plane, i.e. at position (x 0 , y 0 ) = (b, 0).
The remaining injection parameters to investigate are the number of injected deuterium and neon particles, N inj,D and N inj,Ne , and the number of shards, N s,D and N s,Ne . In order to make the most efficient use of the pellet material, the number of shards for a given number of injected particles should be chosen to achieve core penetration without leaving unablated pellet material. To this end, we perform a scan of the  assimilation rate, i.e. the fraction of the pellet material that is deposited in the plasma, as a function of injected quantities and number of shards.
The deuterium assimilation rate as a function of N inj,D and N s,D is shown in figure 2. For a given N inj,D the optimal choice of N s,D is thus close to the 100% assimilation contour. As we will show later, having a good core penetration is essential for a successful disruption mitigation. To ensure a sufficient margin, the number of injected deuterium particles and shards in this paper are chosen along the 97% contour marked by the green line in figure 2.
The relatively low temperature of the diluted plasma when the neon pellets are injected makes their full assimilation difficult. Our studies indicate that the increasing trend in the assimilation rate with increasing number of shards slows down at N s,Ne ∼ 100. Therefore, we fix the number of neon shards to 100.

Representative disruption dynamics
We now consider the spatio-temporal evolution of the most relevant plasma parameters during a representative two-stage injection, using the parameters N inj,D = 2 · 10 24 §, N s,D = 1688, N inj,Ne = 10 23 and N s,Ne = 100. The neon shards are injected at t = 3.4 ms, when the last deuterium shards have just left the plasma centre. We model the diffusive radial electron heat transport due to magnetic perturbations. We switch on a perturbation amplitude of δB/B = 7 · 10 −4 when the neon shards-that create significant pressure variation-enter the plasma, giving a characteristic time scale of the heat transport of about 1 ms. § This quantity roughly corresponds to a 28 mm pellet of the ITER disruption mitigation system.
As the purpose of this study is to investigate trends of disruption mitigation performance measures using a two-stage injection, we consider a simplified case without radial transport of superthermal electrons, reducing the runtime by about an order of magnitude. Such a conservative case can be regarded as an upper limit of the runaway seed generation during the thermal quench, primarily by the hot-tail mechanism. The consideration of such a case is further motivated by the lack of detailed knowledge about the transport of superthermals. The final runaway current has also previously been found to be a logarithmically weak function of the runaway seed during ITER-like conditions [8], and therefore the omission of radial transport of superthermal electrons is not expected to substantially alter the trends found in this study.
The spatio-temporal evolution of the deuterium and neon densities, the electron thermal energy and the current density are shown in figure 3. The two solid green lines mark the trace of the fastest and slowest deuterium shards (for shards with α p = 0), and the dashed green lines mark the corresponding traces for the neon shards. In figure 3a-b, we see a clear increase in the deuterium density and neon density between the corresponding pair of green lines. Note the increase in the neon density saturates quite closely to the first dashed green line, especially in the inner part of the plasma. The neon deposited by the first shards radiatively cools the plasma, impeding the ablation of the later shards.
The sudden deposition of the released thermal energy content during the thermal quench might cause melting of plasma facing components if the heat loads are localised. It is therefore necessary for the disruption mitigation system to ensure that a major part of the thermal energy is lost through radiation. In ITER, the radiated fraction should be larger than 90% of the initial thermal energy content. Looking at figure 3c, we see that the thermal energy density is only slightly affected by the deuterium injection. The plasma is cooled mainly by dilution in this phase, by a factor corresponding to the density increase, resulting in temperatures of a few hundred eV. When the neon shards enter the plasma, the thermal energy density is dissipated-and the temperature dropsover a millisecond time scale or faster in the parts of the plasma that have been reached by the neon shards. The onset of the magnetic perturbations causes some diffusion of the thermal energy density in the parts of the plasma that have not yet been reached by the neon shards. However, in total, diffusive transport dissipates only 7% of the initial thermal energy. Notably, this value is smaller than the 10% limit, and much smaller than the radiative losses that dissipate almost all of the remaining thermal energy. Although it should be emphasised that our model for the thermal energy diffusion is rather crude, assuming here an immediate onset of a prescribed magnetic perturbation that then remains constant in time and space, it allows a simple comparison of the calculated transported losses in different scenarios and can readily be improved for more detailed studies.
The temperature drop to a few eV results in the onset of the current density drop in the part of the plasma which has been reached by the neon shards, as seen in figure 3d. We also see a radial spike in the current profile moving inwards along with the neon shards. This spike is caused by the diffusion of the electric field, induced where the plasma has been cooled, into the hotter region that the neon shards have not yet reached. In this region, the conductivity is still high enough that even a relatively modest increase of the electric field can cause a significant increase in the (ohmic) current density. When the neon shards have reached the core, in most cases, the current starts to decay in all parts of the plasma. In some cases the ohmic heating, amplified by the radial current spike, can cause parts of the plasma to reheat. The re-heating is accompanied by an increase in the conductivity and, as the electric field diffuses into the re-heated regions, the current density can initially increase locally. The decay of the total current is comparatively slow, with a time scale of the order of seconds.

Radiative vs transport losses
As the radiative loss fraction in a disruption is required to be large, we studied the sensitivity of the fraction of the thermal energy lost by transport during twostage injections to the number of injected neon and deuterium particles. In addition to pellet parameters chosen as described in section 3.1, we illustrate the consequences of not achieving core penetration, by including simulations of smaller deuterium pellets which are fully ablated before they reach the core. For these pellets, we chose N s,D = 10. Figure 4 shows the fraction of the initial thermal energy that is lost by transport as a function of the number of injected deuterium and neon particles. We see that the transported fraction can be significantly decreased by increasing the deuterium and neon content.
The dependence on N inj,Ne is weaker than the dependence on N inj,D , mainly due to the decrease in the assimilated fraction as N inj,Ne increases. The neon content primarily increases the radiative losses, which are further enhanced by the increase in the electron density due to the deuterium injection. The deuterium content also limits the transport by lowering the temperature before the onset of the magnetic perturbation. It is noteworthy that simulated transported losses below the limit of 10% of the initial thermal energy, marked by the dashed horizontal line in figure 4, are achievable within a realistic range of injection parameters. Again, further studies with alternative perturbation configurations will be needed to determine the robustness of the observed trend of reduced transported energy losses during two-stage deuterium-neon SPI.
Finally, we note that lack of core penetration does not cause a dramatic difference in the transported fraction. One reason for this is that the neon shards can reach the core even if the deuterium shards did not. Another reason is that the thermal energy in the core has to pass the outer regions of the plasma before it can be lost due to transport, and can therefore be radiated away in the regions that have been reached by the deuterium shards.

Current decay and runaway generation
We now turn our attention to the later stages of the disruption and study the decay of the ohmic current, as well as the runaway generation and dissipation. At this stage, the thermal energy is almost completely dissipated, and the already low plasma temperature evolves slowly as the ohmic heating varies. When the rapid plasma cooling is complete, flux surfaces are expected to re-heal, which we account for by switching off the transport due to magnetic perturbations. The duration of the current quench must be limited to around 50 − 150 ms to limit the forces experienced by the vessel.
We consider separately the non-activated and the nuclear operation phases, with the tritium decay and Compton scattering seed mechanisms being active in the latter. These mechanisms generate a seed runaway current of the order of 0.1−1 A rather independently of the injection parameters. These seeds are sufficiently large to be multiplied to a final runaway current of several MA by the avalanche mechanism, as we will show.
In the relatively cold plasma undergoing a current quench, the radiation transport properties of the plasma can significantly impact the runaway generation. In particular, a preliminary estimate of opacity to Lyman radiation was presented in reference [8], indicating that this effect could lower the runaway currents by up to several MA for large injected deuterium densities.
Not only does the opacity impact the radiative losses, but it also affects the ionization and recombination rates [25]. While the plasma remains essentially transparent at most wavelengths, the opacity is significantly increased at the wavelength corresponding to resonant transitions [26]. This applies particularly to those transitions involving the ground state; as this is the state occupied by most ions and atoms. The estimate shown in Appendix B indicates that the plasma may only be transparent to a few percent of the Lyman radiation resulting from excitations from the ground state, which mostly populates the lower excited states. This impacts the contribution from hydrogen isotopes to the radiation rate L ij . On the other hand, the plasma is estimated to be transparent to the majority of the recombination radiation. A substantial part of this radiation consists of a continuum spectrum resulting from free-bound transitions, together with higher order Lyman lines resulting from the de-excitation of the high excited states thus populated. Opacity to neon radiation is not expected to have a strong impact on disruption dynamics [27]. We therefore consider the limiting cases where the plasma is assumed to be completely transparent or completely opaque to Lyman radiation, whilst remaining completely transparent to radiation from species other than hydrogen. For the completely transparent case, the radiation and ionization/recombination rates are taken from the ADAS database for all species. When the plasma is assumed to be opaque to Lyman radiation, these data are instead taken from the AMJUEL database . Figures 5a-b show the maximum of the runaway current for the different scenarios discussed above, indicating an upper bound of the runaway current that may strike the wall. The thin lines show the nonnuclear case, while the thick lines show results for the nuclear case, i.e. with the tritium decay and Compton scattering seed mechanisms included. Figures 5c-d show the corresponding ohmic current quench times, defined as In figure 5a and c the plasma is assumed to be completely transparent, and in figure 5b and d it is completely opaque to Lyman radiation. The first observation in figure 5a-b concerns the cases with the lowest considered N inj,D and N inj,Ne values (leftmost points of solid curves), where core penetration of deuterium is not achieved. In these cases the radiative cooling is not strong enough to overcome the ohmic heating. This leads to a re-heating of the plasma from the 10 eV range to a few hundred eV, once the transport losses are no longer active. The reheating greatly increases the conductivity, leading to a major reduction of the current decay rate, as well as the induced electric field and the runaway generation rate. In these cases the runaway currents remain relatively small, however the current quench times are unacceptably long (outside the plotted range in panels c-d). The current quench was not completed within the 150 ms of the simulation, but the decay rates indicate a current quench time scale of seconds.
In all other cases shown in figure 5c-d, the current quench times are in the vicinity of the lower acceptable limit of 50 ms, marked by dashed grey lines. Note, however, that in cases with a large runaway conversion, the conversion to runaway current aborts the current quench rather abruptly. The ohmic current quench time calculated here is therefore a lower estimate.
Apart from these, somewhat singular, cases with the lowest injected quantities, the general trend of the maximum runaway current with an increasing N inj,D is either non-monotonic-a decreasing trend turning into an increasing one-in case of a transparent plasma, or a monotonic decrease in case of a plasma opaque to Lyman radiation ¶. The decreasing trend is caused by the following effects. At the postthermal quench temperatures of a few eV, neon is not fully ionized, while deuterium remains practically http://www.eirene.de/html/amjuel.html ¶ Note this the non-monotonic behaviour also appears in a plasma opaque to Lyman radiation, just at higher injected deuterium quantities, outside the range shown here. fully ionized until the temperature drops below ∼ 2 eV. Thus, the injected deuterium contributes with a large increase in the electron density, increasing the critical electric field for runaway generation. Moreover, since a high fraction of bound electrons favours avalanche generation, and that the deuterium content decreases the fraction of bound electrons, we find that avalanche is suppressed when the injected deuterium quantity is increased.
The non-monotonic trend observed in figure 5a is explained as follows. As the ohmic current decays and the ohmic heating decreases, the temperature eventually falls below 2 eV and deuterium starts to recombine. This leads to an increase of the fraction of bound electrons, enhancing the avalanche. At high deuterium densities, the increased radiative losses can cause this to happen, already when there is still a substantial part of the ohmic current left that can be converted to runaways. How much ohmic current remains when the deuterium starts to recombine depends on the radiation transport properties of the plasma. Deuterium recombination also contributes significantly to the radiation, thus opacity to Lyman radiation can reduce the radiative losses considerably. As a consequence, with opacity the ohmic current is smaller when deuterium recombines, hence the lower maximum runaway currents at high deuterium densities in figure 5b. The effect of the radiative properties of the plasma are explored further in the next subsection.
Even though a sufficiently large deuterium injection strongly reduces the hot-tail seed, when tritium decay and Compton scattering is included (thick lines in figure 5a-b) the maximum runaway can not be reduced much below 4 MA for any combination of injected quantities (apart from the case with the unacceptably long t CQ ) even in the presence of opacity effects. This is due to the several orders of magnitude difference between the 0.1 − 1 A seed produced by tritium decay and Compton scattering, and the hottail seed. The finite runaway current depends only logarithmically on the seed, as found in reference [8], due to self-regulating interaction between the runaway current and the electric field. When the runaway current becomes comparable to the remaining ohmic current, the induced electric field decreases, reducing the avalanche growth rate.

Opacity to Lyman radiation
The differences in the dynamics between the transparent and Lyman opaque case are primarily caused by differences in the balance between ohmic heating and radiative losses during the current quench. This balance can be understood by comparing ohmic heating to radiative losses at different temperatures, assuming an  equilibrium distribution of charge states (which may be calculated using equation (7)). The relevant terms, corresponding to the volume averaged final densities of an injection with N D,inj = 2 · 10 24 and N Ne,inj = 10 23 , are plotted in figure 6. The solid black line shows the radiative losses assuming a completely transparent plasma. The green dashed line shows the corresponding radiative losses when the plasma is assumed to be opaque to Lyman radiation. The ohmic heating is calculated for two different current densities: 3.5 MA/m 2 (blue dashed), taken as a representative value for the maximum current density, comparable to the peak value seen in figure 3f, and 0.35 MA/m 2 (blue dotted), representing the phase where the ohmic current has partially decayed. High deuterium density corresponds to strong radiative losses at low temperatures, especially below ∼ 2 eV, where a substantial fraction of the deuterium recombines.
Through recombination, deuterium directly contributes to radiation losses, rather than merely increasing the electron density. In a transparent plasma with j Ohm = 0.35 MA/m 2 , the equilibrium temperature is close to 1 eV, with the ionized fraction of deuterium being only a few percent. This situation favours avalanche, as well as larger induced electric fields (albeit the field decays more rapidly). Moreover, although the ohmic current density is modest in the ∼ 1 eV regions, the conversion from ohmic to runaway current in a given radial location is not limited by the local ohmic current density. The current from other parts of the plasma can diffuse into the cold regions, potentially causing a local increase in the current density. These effects combined give rise to the increase in the maximum runaway current at large deuterium densities seen in figure 5a.
If the plasma is opaque to Lyman radiation, however, the radiative losses at low temperatures are considerably reduced (compare dash-dotted to solid line in figure 6). This is particularly true for the line radiation following excitations from the ground state. The remaining increase in the radiative losses below ∼ 2 eV is instead primarily due to recombination radiation of hydrogen species. The drop to ∼ 1 eV is now postponed until the ohmic current has decreased further than in the transparent case. The postponed temperature drop to ∼ 1 eV reduces the maximum runaway current at large deuterium densities, and shifts the trend of increasing runaway current to even higher deuterium densities.

Discussion
The SPI model used in this work contains a number of simplifications. One such simplification is the use of the NGS ablation model, which is based on a simplified spherical pellet shard geometry and neglecting the details of the ablating electron momentum distribution. In addition, it neglects the electrostatic shielding [28] resulting from the difference in ion and electron mobility and the magnetic shielding [29] due to the deflection of magnetic field lines around the ablation cloud. Nevertheless, the NGS model has been shown to compare reasonably well with experiments, and the effect of the above simplifications have been estimated in the literature to counteract each other [30]. Therefore, while a more advanced ablation model could result in an order unity correction to the deposition profile, the discrepancy compared to the NGS model is expected to be quite modest. A larger correction might be caused by the details of the expansion process of the pellet material between the ablation and final deposition. In the present model, this process has simply been assumed to be instantaneous and local. In reality, however, the homogenisation and equilibration process takes of the order of 1 ms. This time scale is similar to the time it takes the plume of shards to pass a given flux surface. As the newly ablated material is subject to an E × B-drift towards the low field side, the deposited material will cover a wider radial region. (The drift may be less pronounced for neon than for deuterium pellets. The drive of the drift -the excess pressure of the pellet cloud -is directly reduced by the lower temperature of the cloud due to a radiative cooling of neon. In addition, in this two-stage injection scheme the background temperature is lower when the neon shards enter the plasma, with a correspondingly low ablation rate, leading to a lower excess density in the cloud [31].) Apart from the direct impact on the timing and position of the density increase, the details of the expansion process might significantly alter the ablation itself. For pellets with a significant impact on the plasma temperature, the interaction between the plasma and the pellet material is self-regulating; when the pellet material is ablated, the plasma is cooled, at first primarily by dilution, resulting in a slowing down of the ablation. The finite expansion time and drifts alter this self-regulation by delaying the plasma response, and shifting it away from the position of the shards. The ablation might thus be faster, so that more material is deposited earlier along the shard trajectories.
A second area of simplification employed in the present work concerns the geometry and interaction with the structures surrounding the plasma. The geometrical simplifications include the circular plasma cross section, neglect of the toroidicity and the assumption of flux surface homogenised quantities. Relaxing these assumptions would allow for modelling of transient 3D features of the plasma profiles, as well as introducing geometrical order unity corrections to the transport processes involved. For instance, elongating the plasma increases the cross sectional area. This leads to an increase of the time scale for diffusive transport across the plasma cross section, and also decreases the density for a given number of deposited particles. Regarding the surrounding structures, their geometry and conductive properties introduce corrections to the electric field boundary condition. Support for a shaped geometry and a finite wall conductivity are implemented in Dream, and the sensitivity to these features could therefore be studied in the future.
Finally, another simplification in the present model is the prescribed evolution of the magnetic perturbations. The prescribed magnetic perturbation is sufficient to study qualitative trends involving transport due to magnetic perturbations, as done in this work. However, self-consistent and quantitatively accurate simulations would require coupling to a magnetohydrodynamic (MHD) model, such as JOREK [6,32].
The evolution of the magnetic field through the current quench might also be of interest, potentially including deliberately induced perturbations to increase the runaway losses. Even the comparatively smaller magnetic perturbations that may be present during the current quench can have a relatively large impact on the runaway generation and dissipation [33].
This applies especially for cases with an off-axis runaway current profile, as magnetic perturbations induced during the current quench are expected to only partially penetrate into the plasma.
A major drawback with involving three-dimensional MHD modelling is the significant amount of computational resources required. The thermal quench simulations shown in this work take at most a couple of hours on a desktop computer, even with kinetic electrons, while fluid simulations take a couple of minutes. The orders of magnitude lower computational expense substantially increases the feasibility of exploring a wide range of injection and plasma parameters. If the transport processes observed in MHD simulations could be distilled to simplified models suitable for integration in the presented framework, that would allow a major step towards a self-consistent and reliable, as well as computationally efficient, modelling capability for disruption mitigation schemes.

Conclusions
An effective disruption mitigation system limits the exposure of the wall to localized transported heat losses and to the impact of high-current runaway electron beams, as well as avoiding excessive forces on the structure. This work evaluates a two-stage deuterium-neon shattered pellet injection scheme, quantifying these three aspects through the maximum runaway current, the thermal quench timescale, and the transported-to-radiated fraction of thermal energy.
The study is based on simulations of two-stage SPI in ITER-like plasmas using the Dream framework, where an initial dilution cooling achieved by deuterium injection is followed by a radiative thermal quench through neon injection. We find that, when the magnetic field breakup is only triggered by the neon pellet, this injection scheme effectively reduces the contribution of transport to the thermal quench, thereby minimizing localised heat loads. This is due to the lower thermal transport rate in the dilutively pre-cooled plasma.
If the injected neon and deuterium quantities are too low, the current quench time becomes longer than the limit posed by mechanical forces on the wall due to halo currents. The long ohmic current quench time is caused by an incomplete cooling in parts of the plasma, and the associated high conductivity. For large injected densities, however, the simulated current quench times were found to be close to the lower acceptable limit, with a rather weak dependence on the injected densities.
We found that the two-stage SPI scheme can reduce the hot-tail runaway seed generation by several orders of magnitude. This reduction is explained by the thermalization of the electron distribution between the injections.
The maximum runaway current exhibits a nonmonotonic dependence on the injected deuterium density: At low deuterium densities, an increased deuterium density reduces the avalanche, while at large deuterium densities, the avalanche is enhanced by the higher electric fields and substantial recombination, resulting from the strong radiative cooling. An optimum is found at around an injected deuterium quantity of 1 − 2 · 10 24 atoms, depending on the opacity of the plasma to Lyman radiation. Subsequent neon injection in the range of 10 22 − 10 24 atoms gives current quench times and transported fractions of the thermal energy within/close to the respective ranges of acceptable values.
In nuclear operation no combination of injected quantities could reduce the maximum runaway current much below 4 MA, which is an alarming result. The runaway current carried by the beam upon wall impact might be significantly affected in the presence of naturally occurring or externally applied magnetic perturbations remaining after the thermal quench [32], which was not considered in this work.
Although the model used in this work comprises an integrated framework accounting for many of the relevant aspects of disruption mitigation by SPI, the various components are treated with different levels of sophistication. For instance, heat and particle transport processes are accounted for in a simplified manner, while runaway electron dynamics is modelled in a comprehensive and state-of-the-art fashion. The benefits of a two-stage SPI scheme indicated by our results motivate further quantitatively accurate studies. No 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

Appendix A. Overview of Dream
Dream solves for the flux-surface averaged electron and ion densities, temperatures, as well as the parallel electric field and current. The electron dynamics is resolved fully kinetically or through various fluid models. Going beyond the status of the Dream framework as represented in reference [11], here we overview the aspects of the code relevant for SPI modelling.

Particle and energy balance
In the simulations presented here the electron distribution is resolved kinetically for the thermal and superthermal electron populations, while the highly energetic runaway electrons-above a momentum p RE which we here set to 3m e c-are treated as a fluid. Even when the thermal part of the electrons is treated kinetically, fluid energy conservation and quasineutrality equations are also being solved for a Maxwellian bulk electron species to account for radial transport and atomic physics processes. The kinetic collision operator is linearized around this fluid bulk, which enforces that the kinetically resolved thermal electron population remains close to the fluid one, and electron sources on the kinetic grid make sure that the fluid and kinetic particle numbers are consistent.
The evolution of the energy density W M = 3n M T M /2 of the Maxwellian bulk electrons is governed by The conductivity σ M appearing in the ohmic heating term is calculated using the relativistic expression derived in reference [34] σ M =σ whereσ(T M , n M , Z eff ) is calculated by interpolation of the values tabulated in reference [34]. Here we have also introduced the effective charge Z eff = ij Z 2 ij n ij /n free , the dielectric constant 0 , the elementary charge e, the thermal Coulomb logarithm ln Λ 0 = 14.9 − 0.5 ln (n M /10 20 ) + ln (T M /10 3 ) [35], and the electric field parallel to the magnetic field lines E || . The second term accounts for the energy required to ionize the ablated pellet material and reach the equilibrium charge state distribution is the total energy required to ionize an atom of species j from neutral to charge state i, and the ionization energies E ioniz ij are taken from the NIST database + . The assumption of instantaneous equilibration and homogenization over the flux surface, of the ablated material, means that the thermal energy absorbed by the shielding cloud is immediately returned to the background plasma. We therefore do not need any further energy loss terms directly associated with the pellet ablation (assuming the shielding cloud to be optically thick, thus neglecting radiative losses from it).
The third term in equation (A.1) describes electron heat diffusion. Here we employ a Rechester-Rosenbluth-type [36] diffusion coefficient D = πqv || R 0 (δB/B) 2 , for particles with a parallel streaming speed of v , where the relative amplitude and the parallel correlation length of magnetic perturbations are δB/B and πqR 0 , respectively, with q ≈ 1. Integrating over the Maxwellian electron bulk yields where v T = 2T M /m e is the electron thermal speed.
Here, whenever we switch on transport, we prescribe δB/B = 7 · 10 −4 . With D W evaluated at the initial temperature, this results in a transport loss time scale of a Bessel mode-like decay a 2 /(D W x 2 1 ), with x 1 ≈ 2.4, + https://physics.nist.gov/PhysRefData/ASD/ionEnergy. html comparable to the expected thermal quench time in ITER [37]. The fourth term in equation (A.1) corresponds to the line radiation and ionization losses. The line radiation rates L ij (T M , n M ) are taken from the OpenADAS database [24]. The charge state distribution is calculated by the time dependent rate equations Thus, the total evolution of the ion charge state densities is given by The electron density is determined from the quasineutrality condition. The fifth term in equation (A.1) describes the collisional energy transfer from hot electrons to the Maxwellian electrons * .
Here ∆Ė ee = 4πn M r 2 0 ln Λ ee m e c 4 /v, where r 0 = e 2 /(4πε 0 m e c 2 ) is the classical electron radius and ln Λ ee is the energy dependent Coulomb logarithm for electron-electron collisions, given in reference [38]. The last term in equation (A.1) accounts for the bremsstrahlung losses.

Electric field evolution
In the cylindrical limit, Faraday's law combined with Ampère's law yield where j || is the parallel current density, given by the sum of the ohmic current density, and the hot electron and runaway current densities. The ohmic current density is calculated as j Ohm = E || (σ M − σ fp ) + p hot 0 ev || f dp, where σ M is the conductivity given in equation (A.2). The first term corrects for the part of the conductivity not captured by the conductivity σ fp resulting from the test-particle Fokker-Planck collision operator used in the kinetic equation. An expression for σ fp was determined by running numerous DREAM simulations with fixed parameters until the kinetically captured contribution to the ohmic current was equilibrated, and then calculating σ fp by dividing the ohmic current thus obtained by the fixed value for E || . This data was used to fit an expression for σ fp according to , (A.8) * As this depends on the energy distribution, which is not resolved for the runaway population; here, only the hot population is accounted for in this interaction term.
with c 1 = −1.406 and c 2 = 1.888. The hot electron current density is j hot = p hot <p<pRE ev || f dp, and the runaway current density is j RE = ecn RE . The boundary condition for equation (A.7) at r = a is obtained by assuming the plasma to be surrounded by a perfectly conducting wall at r = b > a, where the induced electric field is set to 0. Matching the solution for r < a to the vacuum solution for a < r < b gives E || (a) = a ln (a/b)∂E || /∂r| r=a .

Numerical resolution
The kinetically resolved electron momentum space, that covers momenta up to the boundary of the numerical runaway region p RE = 3m e c, is divided into two regions: Below p = 0.07m e c, covering the thermal bulk, the momenta are resolved with 70 uniform grid cells. Above this, and below p RE , momenta are resolved with 50 uniform grid cells. The pitch angle cosine ξ is resolved by 5 grid cells uniformly spaced between −1 and 1 in the entire kinetically resolved region. The radial dimension is covered uniformly by 11 grid cells. The temporal resolution varies between 1 − 10 µs, using values from the lower end of this range in the thermal quench phase. The relative tolerance of unknowns, when advancing the system one time step by a nonlinear iterative solver, is set to 10 −6 .

Appendix B. Opacity for deuterium radiation
Here we estimate the fraction of the deuterium line radiation which is trapped due to opacity of the plasma in the vicinity of the Lyman lines. We consider a plane, partially ionized, plasma slab, and follow a model described in references [26] and [27]. The fraction of trapped radiation is determined mainly by two quantities: the optical thickness of the plasma to the line radiation, and the rate of collisional quenching of excited states. If the rate of collisional quenching is low, even an optically thick plasma can be affected by strong radiative energy losses, as the absorption of any photons would quickly be followed by a new photon emission, leading to an efficient radiative transport. The optical thickness of a plasma slab of thickness h is given by τ = k l ij h, with k l ij being the inverse mean free path of the photon emitted by deexcitation from the l th excited state to the ground state of charge state i of ion species j (transitions between excited states are neglected). We will limit ourselves to the first nine excited states due to availability of the necessary data (providing an estimated accuracy of ∼ 1%). With the natural broadening of the line profile γ and the net external broadening Γ, the inverse photon mean free path is Here, n ij is the density of charge state i of ion species j and λ l ij is the wavelength of a photon emitted by deexcitation from the l th excited state to the ground state of charge state i of ion species j. For a Doppler broadened line, we have Γ/γ = 1.11 · 10 10 E 0l ij m p T I /m I ν l ij , where m p is the proton mass, m I is the ion (or atom) mass, T I is the ion temperature, E 0l ij is the energy difference between the ground state and the l th excited state of charge state i of ion species j and ν l ij is the corresponding natural decay rate. Here, all temperatures and energies are given in eV, and transition rates in s −1 .
The effect of collisional quenching of excited states is determined by the ratio of the collisional quenching probability and the decay probability, where E 1 (x) = ∞ x dt exp (−t)/t is the integral exponent. The fraction of radiation escaping from a plane ion slab for a line emitted by deexcitation of the l th excited state from charge state i of ion species j is then given by B l ij = W l ij /(β l ij + 1), where W l ij = (1 + β l ij )P a /(β l ij + P a ), and P a = [1 + τ π ln (τ + 1)] −1 is the probability for a photon to travel a distance h without being absorbed (the "1" in the denominator is included ad-hoc in the expression valid in the optically thick limit, to make sure that P a → 1 as it should for thin plasma slabs).
To find the fraction of the total radiation that escapes the plasma, we also need to know the relative intensity of the different lines. The intensity distribution over the different lines is very different for excited states populated by excitations from lower states (dominantly the ground state) compared to excited states populated by recombination. The reason for this is that excitations from lower states primarily populate the lower excited states, while recombination populates higher excited states. In addition, part of the potential energy change is radiated away during a freebound transition, resulting in a continuous spectrum to which the plasma is essentially transparent.
For radiation following excitation from lower states, we calculate the relative line intensities based on data from reference [39], according to L l ij ∝ E 0l ij φ l ij , where φ l ij = r l ijn Saha (l)ν l ij = r l ij (l + 1) 2 exp (E l∞ ij /T e )ν l ij is proportional to the relative occupationn Saha of the l th excited state at Saha equilibrium, and E l∞ ij is the ionization energy from the l th excited state. The coefficients r l ij (tabulated in [39]) describe the deviation from the Saha equilibrium prevailing in a dense enough plasma to reach local thermodynamic equilibrium. The transition rates ν l ij are tabulated in [40]. The fraction of line radiation following excitations from the ground state escaping the plasma can now be expressed as where the second equality holds for hydrogen isotopes where there is only one radiating charge state. When the excited states are populated by recombination, energy will also be radiated away during the free-bound transition, in addition to the subsequent deexcitation through bound-bound transitions. We consider two types of recombination event: one with transitions directly from a free state to the bound ground state, and one consisting of a free-bound transition to an excited bound state, followed by a direct transition to the ground state. If the recombination is radiative, a transition from state l (which can be the ground state, l = 0) will be preceded by a transition from the continuum to state l, while releasing an average photon of energy E l ij,rec ≈ E ∞0 ij + T e − E 0l ij . The free-bound transitions give rise to a continuous spectrum, to which the plasma is essentially transparent. The resulting escape factor for recombination radiation becomes f esc = ij n ij l (B l ij L l ij + L l ij,rec ) ij n ij l (L l ij + L l ij,rec ) = l (B l 0D E 0l 0D + E l 0D,rec )φ l 0D l (E 0∞ 0D + T e )φ l 0D .

(B.4)
As the coefficients in [39] only apply to recombination events involving an excited bound state, we replace φ l 0 in equation (B.4) with the transition rates from [41], which also include radiative recombination directly to the ground state. The escape factor as a function of the slab thickness h is shown in figure B1, for h up to the minor radius of an ITER-like plasma a = 2 m. Results are shown for both radiation following excitations from the ground state (solid), from equation (B.3), and for recombination radiation (dotted), from equation (B.4). The plasma parameters considered, n e = 10 20 m −3 and T e = 1.38 eV, and neutral deuterium density n D = 4 · 10 21 m −3 , are chosen to match the ones for which there are tabulated data in reference [39]. One can see that the line radiation following excitations is below 1% already for h = 0.5 m. For the recombination radiation, on the other hand, the escaping fraction remains above 60% for the h-range shown.  Figure B1. Escaping fraction of deuterium radiation as a function of slab thickness, for a plasma with electron density ne = 10 20 m −3 , neutral deuterium density n D = 4 · 10 21 m −3 , and temperature Te = 1.38 eV. Black solid: Radiation following excitation from the ground state. Dotted red: recombination radiation.