# Collisional-radiative simulation of impurity assimilation, radiative collapse and MHD dynamics after ITER shattered pellet injection

D. Hu, E. Nardon, F.J. Artola, M. Lehnen, D. Bonfiglio, M. Hoelzl, G.T.A.

Huijsmans, S.-J. Lee

► To cite this version:

D. Hu, E. Nardon, F.J. Artola, M. Lehnen, D. Bonfiglio, et al. Collisional-radiative simulation of impurity assimilation, radiative collapse and MHD dynamics after ITER shattered pellet injection. Nuclear Fusion, 2023, 63 (6), pp.066008. 10.1088/1741-4326/acc8e9. cea-04707701

HAL Id: cea-04707701

https://cea.hal.science/cea-04707701v1

Submitted on 24 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Collisional-radiative simulation of impurity
assimilation, radiative collapse and MHD dynamics
after ITER Shattered Pellet Injection

D. Hu<fnref n="1" />, E. Nardon<fnref n="2" />, F. J. Artola<fnref n="3" />, M. Lehnen<fnref n="3" />, D. Bonfiglio<fnref n="4" />,
M. Hoelzl<fnref n="5" />, G.T.A. Huijsmans<fnref n="2" />,<fnref n="6" />, S.-J. Lee<fnref n="7" /> & JOREK Team‡

E-mail: hudi2@buaa.edu.cn

15 March 2023

**Abstract.** Recent studies suggest significant time delay between the Shattered Pellet Injection (SPI) fragment arrival and the temperature radiative collapse could exist in ITER, depending on the impurity assimilation and the plasma thermal reservoir. Hence in some cases the fragments could reach the core even before the edge radiative collapse occurs and triggers strong stochastic transport. This could be beneficial for heat load mitigation and hot-tail runaway electron suppression. To investigate the expected assimilation and radiation, thus the MHD response after SPIs in 3D, we carry out simulations of collisional-radiative impurity mixed SPIs into ITER L-mode equilibrium. Localized cooling around the fragments is found to cause current perturbations which destabilize MHD modes. Meanwhile, slower injections are found to result in stronger and more complete radiative collapse, thus stronger MHD amplitude. Due to the $q = 1$ surface enclosing a significant volume, the 1/1 resistive kink mode is shown to couple with outer modes to bring global stochasticity and convective core density mixing, although a transport barrier outside of the $q = 1$ surface prevents immediate temperature relaxation over the whole plasma. The impact of various physical assumptions and numerical treatments, such as the use of the flux-averaged ambient plasma parameters for ablation calculation, the exclusion of the magnetic constraining effect in ablation, the localization of the density source and the use of constant parallel thermal conduction instead of the Braginskii one and different injection velocities are also investigated. In general, stronger and more localized ablation results in stronger radiation, faster radiative collapse and a more violent MHD response, while the assimilation changes little due to a self-regulation effect.

<footnote n="1">¹Beihang University, No. 37 Xueyuan Road, Haidian District, 100191 Beijing, China.</footnote>

<footnote n="2">²CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France</footnote>

<footnote n="3">³ITER Organization, Route de Vinon sur Verdon, CS 90 046,13067 Saint Paul-lez-Durance, Cedex, France.</footnote>

<footnote n="4">⁴Consorzio RFX-CNR, ENEA, INFN, Università di Padova, Acciaierie Venete SpA.
I-35127 Padova, Italy.</footnote>

<footnote n="5">⁵Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching b. M., Germany</footnote>

<footnote n="6">⁶Eindhoven University of Technology, De Rondom 70 5612 AP Eindhoven, the Netherlands.</footnote>

<footnote n="7">⁷Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, South Korea.</footnote>
# 1. Introduction

The inward-propagating radiative collapse of the electron temperature profile is an important and characteristic behaviour after impurity mixed Shattered Pellet Injection (SPI) found both in experiments [1–4] and simulations [5–7]. Such drastic temperature collapse would result in significant resistive current density redistribution and consequently strong current driven magneto-hydrodynamic (MHD) instabilities, leading to the onset of a mitigated Thermal Quench (TQ) [6]. Such current redistribution could either be axisymmetric [8, 9], or helical [10–13], or a combination of both. In present devices, the electron temperature collapse is usually so complete and fast, that remarkable current redistribution closely follows the impurity-mixed fragments, making the pre-TQ dynamics of even very small mixture ratio impurity/deuterium SPI different from that of a pure deuterium one [2, 4].

However, recent 1.5D and 2D simulations have shown that a significant time delay could exist between the arrival of mixed impurity fragments and the occurrence of edge temperature radiative collapse for future large devices such as ITER [14, 15]. The exact delay time depends on various parameters such as the injection velocity, the injection assimilation and the plasma thermal reservoir. For example, it is found in those numerical investigations that, with 500m/s injection velocity, 5% neon mixture ratio and reasonable fragmentation size, it takes about 10ms for the radiative collapse to occur around the $q = 2$ surface for ITER L-mode target equilibria, longer than the time for fragments to reach the magnetic axis [14]. An electron temperature plateau on $O(100eV)$, due to dilution cooling, is also found in the wake of the fragments before the radiative collapse ultimately occurs [14]. On the other hand, with experimental SPI parameters, the radiative collapse around the $q = 2$ surface is found to occur in about 1 millisecond for DIII-D plasmas, in good agreement with experimental observation [15]. No electron temperature plateau is found before the radiative collapse in this case [15]. Such difference in the pre-TQ dynamics suggests some “size effects” may play important roles in determining the pre-TQ and the ensuing early TQ dynamics.

Although already informative, the above investigations did not include 3D effects such as the coupling with MHD activities which incurs significant density mixing and temperature relaxation [6, 13], as well as localization effect of the impurity radiation sink before the ablation cloud finally relaxes in finite time [16, 17], which should be taken into account for a more comprehensive albeit more expensive study. To take these factors into account and to provide more insights into the pre-TQ and early TQ cooling dynamics in ITER, we hereby use the 3D nonlinear MHD code JOREK [18] to simulate the cooling, MHD response and the ensuing density transport after neon mixed SPI into ITER L-mode equilibria. The focus of this study is on the pre-TQ and early TQ cooling dynamics and its impact on the MHD response, thus we do not delve too much into other important issues such as the radiative fraction and the radiation asymmetry etc., which will be explored in future works. The reduced MHD version of the code will be used in this study.
The rest of the paper is arranged as the following: In Section 2, our system of interest is introduced. We first introduce our governing reduced MHD equations and noteworthy numerical treatments. We then proceed to introduce the target equilibria, simulation setup as well as the SPI configurations and parameters. Then in Section 3, the characteristic temperature cooling, MHD dynamics and the ensuing density transport after a neon SPI injection are shown, and the transient electron temperature behaviour referred to as the “electron temperature hole” is described and their impact on the MHD response is discussed. The impact of injection velocities on the pre-TQ cooling dynamics is shown. The transport barrier close to the $q = 1$ surface is also discussed. We then further examine the impact of various assumptions and numerical treatment in Section 4, with focus on the difference in impurity assimilation, occurrence of temperature collapse and the radiated energy. Finally, conclusions and discussions on the expected pre-TQ and early TQ dynamics after impurity SPI in ITER are presented in Section 5.

## 2. The system of interest

The system of interest for this study is introduced in this section. In Section 2.1, we briefly introduce our governing equations and basic assumptions used in our JOREK simulations, and our ITER L-mode target equilibrium is introduced in Section 2.2. The SPI configuration and parameters are shown in Section 2.3. Last, in Section 2.4, we present a discussion on the so-called “size effects” and its expected impact on the pre-TQ dynamics.

### 2.1. The governing equations and assumptions

The reduced MHD version of the JOREK 3D nonlinear code [18] is used to carry out simulations in the following sections. We use the two-temperature version of the collisional-radiative non-Coronal-Equilibrium impurity model detailed in Ref. [6] and Ref. [19] where ionization and excitation are dominated by collision while recombination and line-radiation are dominated by radiative process [20]. We briefly write down the governing equations in the tokamak coordinates $(R, Z, \phi)$ for the completeness of the paper, without going into great details:

$$ \mathbf{B} = F_0 \nabla \phi + \nabla \psi \times \nabla \phi, \qquad (1) $$

$$ \mathbf{v} = v_{\parallel} \mathbf{B} - R^2 \nabla u \times \nabla \phi, \qquad (2) $$

$$ j = \Delta^* \psi \equiv R^2 \nabla \cdot \left( \frac{1}{R^2} \nabla_{pol} \psi \right), \quad j_{\phi} = -j/R, \qquad (3) $$

$$ \omega = \frac{1}{R} \frac{\partial}{\partial R} \left( R \frac{\partial u}{\partial R} \right) + \frac{\partial^2 u}{\partial Z^2}, \qquad (4) $$

$$ \frac{\partial \psi}{\partial t} = \eta(T_e) \Delta^* \psi - R \{u, \psi\} - F_0 \frac{\partial u}{\partial \phi}, \qquad (5) $$
$$ \frac{\partial \rho}{\partial t} = - \nabla \cdot (\rho \mathbf{v}) + \nabla \cdot (D \nabla \rho) + S_{bg} + S_{imp}, \qquad (6) $$

$$ \frac{\partial \rho_{imp}}{\partial t} = - \nabla \cdot (\rho_{imp} \mathbf{v}) + \nabla \cdot (D \nabla \rho_{imp}) + S_{imp}. \qquad (7) $$

$$ R \nabla \cdot \left[ R^2 \frac{\partial}{\partial t} (\rho \nabla_{pol} u) \right] = \frac{1}{2} \left\{ \begin{aligned} & R^2 |\nabla_{pol} u|^2, R^2 \rho \end{aligned} \right\} + \left\{ \begin{aligned} & R^4 \rho \omega, u \end{aligned} \right\} \\ & - R \nabla \cdot \left[ R^2 \nabla_{pol} u \nabla \cdot (\rho \mathbf{v}) \right] + \{\psi, j\} \\ & - \frac{F_0}{R} \frac{\partial j}{\partial \phi} + \{P, R^2\} + R^3 \mu_{\perp} (T_e) \nabla_{pol}^2 \omega, \qquad (8) $$

$$ B^2 \frac{\partial}{\partial t} (\rho v_{\|}) = - \frac{1}{2} \rho \frac{F_0}{R^2} \frac{\partial}{\partial \phi} (v_{\|} B)^2 - \frac{\rho}{2R} \{B^2 v_{\|}^2, \psi\} - \frac{F_0}{R^2} \frac{\partial P}{\partial \phi} \\ + \frac{1}{R} \{\psi, P\} - B^2 \nabla \cdot (\rho \mathbf{v}) v_{\|} + B^2 \mu_{\|} \nabla_{pol}^2 v_{\|}. \qquad (9) $$

$$ \frac{\partial}{\partial t} P_i = - \mathbf{v} \cdot \nabla P_i - \gamma P_i \nabla \cdot \mathbf{v} + \nabla \cdot (\kappa_{\perp} \nabla_{\perp} T_i + \kappa_{i, \|} \nabla_{\|} T_i) \\ + \frac{\gamma - 1}{2} \mathbf{v} \cdot \mathbf{v} (S_{bg} + S_{imp}) + (\gamma - 1) \mu_{\|} [\nabla_{pol} (v_{\|} B)]^2 \\ + (n_{bg} + n_{imp}) (\partial_t T_i)_{c,e}, \qquad (10) $$

$$ \frac{\partial}{\partial t} P_e = - \mathbf{v} \cdot \nabla P_e - \gamma P_e \nabla \cdot \mathbf{v} + \nabla \cdot (\kappa_{\perp} \nabla_{\perp} T_e + \kappa_{e, \|} \nabla_{\|} T_e) \\ + \frac{\gamma - 1}{R^2} \eta (T_e) j^2 - (\gamma - 1) (P_{rad} + P_{ion}) + n_e (\partial_t T_e)_{c,i}. \qquad (11) $$

In the above equations, the Poisson bracket $\{f,g\} \equiv R(\nabla f \times \nabla g) \cdot \nabla \phi$ and $\nabla_{pol}$ represents the gradient within the poloidal plane. Eq. (1) and Eq. (2) give the magnetic and velocity fields. Eq. (5) is the induction equation, and $\eta$ is the local Spitzer resistivity [22]. Note that the effective charge contribution [23] is included

$$ \eta = Z_{\text{eff}} \frac{\eta_0}{\max(T_e, T_{\text{thres}})^{3/2}} \times \frac{1 + 1.198 Z_{\text{eff}} + 0.222 Z_{\text{eff}}^2}{1 + 2.966 Z_{\text{eff}} + 0.753 Z_{\text{eff}}^2}. \qquad (12) $$

Here the threshold temperature $T_{thres} = 1eV$, below which the resistivity no longer varies as a function of the temperature.

Further, Eq. 6 and Eq. 7 are the total mass continuity equation and the impurity mass continuity equation which include all impurity charge states from neutrals to fully ionized ions. The background species (hydrogen in this particular study) and impurity ablative density sources $S_{bg}$ and $S_{imp}$ are given by a Gaussian shape distribution around each fragment with position $(R_f, Z_f, \phi_f)$:

$$ S_n \propto \exp \left( - \frac{(R - R_f)^2 + (Z - Z_f)^2}{\Delta r_{NG}^2} \right) \times \exp \left( - \left( \frac{\phi - \phi_f}{\Delta \phi_{NG}} \right)^2 \right). \qquad (13) $$

In this paper we choose $\Delta r_{NG} = 0.08m$. Although the individual density source width from this choice may be under-resolved in the poloidal direction, the total density source
consists of an envelop of all fragment contribution from the plume hence is on the order of $0.5m$, thus the ablation density source is generally well resolved. This also makes the result insensitive to the exact choice of $\Delta r_{NG}$ so long as it is not larger than the plume width. Our toroidal resolution limit prevents us from using realistic $\Delta\phi_{NG}$, as that requires tracking a high amount of toroidal harmonics in our Fourier decomposition and makes the simulations extremely expensive. In order to assess the sensitivity of the dynamics to $\Delta\phi_{NG}$, we test two values: $\Delta\phi_{NG} = 0.3$ and $0.5$ rad. The impact of such artificial toroidal elongation will be discussed in Section 4.2. The total mass ablation rate for each neon-hydrogen mixed fragment is determined by the Neutral Gas Shielding (NGS) ablation model [24]. The NGS model estimates the ablation rate needed to maintain a certain line-integrated density in the neutral cloud around the pellet, so that the parallel heat flux can barely be completely shielded on the pellet surface [25]. This version of NGS model we use also takes into account of the areal deposition of the density source, the magnetic reduction of the ablation rate as well as capable of describing the ablation rate of mixed pellet. The ablation rate for each fragment is explicitly:

$$G[g/s] = \lambda(X) \zeta_B \left(\frac{T_e[eV]}{2000}\right)^{5/3} \left(\frac{r_p[cm]}{0.2}\right)^{4/3} n_e [10^{14}cm^{-3}]^{1/3}, \quad (14)$$

in which the parameter $\lambda$ is a function of the mixture ratio $X \equiv \frac{N_{D2}}{N_{Ne}+N_{D2}}$:

$$\lambda(X) \equiv 27.08 + \tan(1.49X), \qquad (15)$$

$r_p$ is the fragment radius and $\zeta_B \approx (\frac{2}{\max(B[T],2)})^{0.843}$ is a fitted parameter to represent the reduction of ablation in the strong magnetic field regime due to the smaller perpendicular expansion of the plasmoid [24]. The density source for each species is then obtained by considering the pellet mixture ratio. Note that we can choose either to use the local fluid electron temperature $T_e$ and density $n_e$ in calculating the ablation rate, or use the flux-surface-averaged values instead to take into account the long mean-free-path of the hot electrons “seen” by the pellet fragment [25, 26]. For simplicity, we hereby call the two choices the “local ablation” and “non-local ablation” scheme respectively. The physical incentive of such “non-local ablation” scheme is discussed at the end of Section 3.2 and its impact is discussed later on in Section 4.1.

Moreover, the perpendicular and the parallel momentum equations are shown in Eq. (8) - Eq. (9), and we have assumed that all species and all charge states share the same velocity due to frequent collision or charge exchange. At last, Eq. (10) and Eq. (11) are the ion and electron pressure equations. The Braginskii parallel heat conduction [27] is used for each species respectively. With the initial density, the perpendicular particle diffusion is set to be constant at $6m^2/s$ and the perpendicular heat conduction at $10m^2/s$. The perpendicular viscosity is set to $30m^2/s$ and the parallel one to $150m^2/s$. The collisional thermalization between the ions and the electrons is represented by $(\partial_t T_i)_{c,e}$ and $(\partial_t T_e)_{c,i}$ [6]:

$$(\partial_t T_e)_i = (\nu_{e/imp} + \nu_{e/bg}) (T_i - T_e), \quad (\partial_t T_i)_e = -\frac{n_e}{n_{imp} + n_{bg}} (\partial_t T_e)_i,$$
where $\nu_{e/imp}$ and $\nu_{e/bg}$ are the electron-ion collision rate for the impurity and background species respectively. The ionization power density and the radiation power density are $P_{ion}$ and $P_{rad}$ respectively. Note that we assume all the background species neutrals are instantaneously ionized while the impurity ionization, recombination and radiation are described by the collisional-radiative probability matrices from the open ADAS data base [28]. The charge state evolution is then determined by [19]:

$$ \frac{df_i}{dt} = n_e [S_{i-1}(T_e)f_{i-1} - (S_i(T_e)f_i + \alpha_i(T_e))f_i + \alpha_{i+1}f_{i+1}]. \quad (16) $$

Here $f_i$ is the fractional abundance, $S_i$ and $\alpha_i$ are the corresponding ionization and recombination probabilities and $n_e$ is the electron density. For a given charge state distribution, the radiative power is then [21]

$$ \begin{aligned} P_{rad} &= n_e n_{imp} L_{rad}(n_e, T_e) \\ &= n_e \sum_{i=0}^{Z} n_{imp}^{i+} (L_{L}^{i+}(n_e, T_e) + L_{R}^{i+}(n_e, T_e) + L_{B}^{i+}(n_e, T_e)). \end{aligned} \quad (17) $$

Here $n_{imp}$ is the impurity density including all charge states while $n_{imp}^{i+}$ is that of each individual state, $L_{L}^{i+}$, $L_{R}^{i+}$ and $L_{B}^{i+}$ are the corresponding radiation function for the line radiation, the recombination radiation and the bremsstrahlung radiation. It should be noted that all impurity related values such as the effective charge and the impurity ionization and radiation powers are obtained using the projection of the marker-particle moment as is described in Ref. [6].

## 2.2. The target ITER L-mode equilibrium and the simulation setup

An ITER L-mode equilibrium is chosen as the default target plasma [29]. The toroidal magnetic field is $B_T = 5.3T$ and the total plasma current is $I_p = 15MA$, the central pressure is $P_0 = 8.654 \times 10^4 Pa$, the total initial thermal energy is $E_{th}(0) = 36MJ$. The initial equilibrium profiles, including the electron temperature, density, the toroidal current density and the safety factor $q$ profile, are shown in Fig. 1. This equilibrium exhibits a slightly reversed shear near the magnetic axis, and a significant portion of the core plasma (up to normalized flux $\bar{\psi} \simeq 0.3$) is within the $q = 1$ surface, thus we expect a strong 1/1 resistive kink response even if our reduced MHD description tends to underestimate its instability [18, 30]. Without SPI, our numerical investigation does not see significant 1/1 kink growth within our timescale of interest. For simplicity, we set the simulation boundary at the approximate position of the thin wall equivalent of the ITER blanket modules [31]. The grid is shown in Fig. 2. In this study, we trace 12 toroidal harmonics from $n = 0$ to $n = 11$. With this range of toroidal harmonics, we are able to keep Fourier components up to the first order in strength for the $\Delta\phi_{NG} = 0.3$ cases and up to the third order in strength for the $\Delta\phi_{NG} = 0.5$ cases. Ideal boundary condition is assumed in this study. We call this equilibrium the “standard equilibrium”.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The initial equilibrium. $\Psi_N$ is the normalized poloidal flux.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The simulation grid. The boundary is assumed to be ideally conducting.</figcaption></figure>
## 2.3. The SPI configuration

The SPI is carried out from the outer equatorial port along the major radial direction. The vertex angle of the injection cone is 20 degrees, and the velocity spread is ±40%. We assume the fragment velocity distribution function is flat for simplicity, although in reality the velocity distribution has a correlation with the fragment size [32]. This means that we would expect finer fragments with higher injection velocity in reality, enhancing the surface-to-volume ratio for these cases which would result in stronger ablation than predicted here. We do not expect this to strongly alter the dynamics we investigate here, since one could usually modify the shattering angle to partly cancel this effect. We have two choices of fragment plume reference velocity in this study, the “slow case” with $v_{p,0} = 200m/s$ and the “fast case” with $v_{p,0} = 500m/s$. The fragment size distribution follows that of the Statistical Fragmentation model [33, 34]:

$$P(r_p) = \frac{r_p K_0(\kappa_p r_p)}{I}, \quad I = \int_0^\infty r_p K_0(\kappa_p r_p) dr_p = \kappa_p^{-2}. \qquad (18)$$

Here $r_p$ is the radius of the fragments and $P(r_p)$ is the probability distribution of fragment size, $K_0$ is the modified Bessel function of the second kind, $\kappa_p$ is the inverse characteristic fragment radius, which in principle should scale as $L^{-1}$ where L is the length of the un-shattered pellet, for fixed shattering angle, mixture ratio and normal velocity against the shattering plate [34]. This suggests that the characteristic fragment size goes up linearly with the un-shattered pellet size. Integrating Eq. 18 over $r_p$, we get:

$$N_p n_p 6\pi^2 \kappa_p^{-3} = N_0, \qquad (19)$$

where $N_0$ is the total amount of atoms turned into solid fragments, $n_p$ is the atom number density of the solid fragments and $N_p$ is the total number of fragments.

In this study, we consider the injection of $5 \times 10^{22}$ neon atoms and $1.8 \times 10^{24}$ deuterium atoms. Without going into details of shattering angle and normal injection velocity, we assume all the injected mass is converted into 300 solid fragments without gas or droplet generation in this study. The corresponding characteristic fragment size is $\kappa_p^{-1} \simeq 1.43mm$. Note that here we do not consider the correlation between the shattering fineness and the injection velocity as well as the mixture ratio [32, 34] for the sake of simplicity. We also do not consider the velocity change of the fragment due to its interaction with the target plasma such as the rocket effect or the plasma drag force in this study.

## 2.4. Discussion of the “size effects”

Here we discuss the so-called “size effects” regarding the time it takes for a given “layer” of flux surfaces to cool down significantly after the fragmetns first arrive on these surfaces. Assuming the same thermal energy density, large devices such as ITER naturally contain more thermal energy. Taking the aspect ratio as a near constant, for a given “layer” of flux surfaces covering the same portion of normalized minor radius, the
thermal energy contained within goes up with $a^3$ where $a$ is the minor radius. On the other hand, the total injection amount also often goes up with the volume $N_0 \propto V_{plasma}$, so that the thermal energy per injected atom should, naively, be constant. Thus one would crudely anticipate the timescale of the dilution and radiation cooling after SPIs in larger devices not to change much compared with existing machines due to this cancellation.

However, as we mentioned in the last subsection, the characteristic fragment radius also goes up with the pellet length: $\kappa_p^{-1} \propto L$. Combined with Eq. 19, this means that, assuming the same pellet shape, shattering angle, injection velocity and mixture ratio, the total number of solid fragments $N_{frag}$ should remain more or less the same even with changing un-shattered pellet size. Meanwhile, with the same $T_e$ and $n_e$, the ablation rate goes up with $r_p^{4/3}$, so that the total ablation rate from all fragments should have the following relationship with the injection amount: $G_{tot} \propto B^{-0.842} \kappa_p^{-4/3} \propto B^{-0.842} L^{4/3} \propto N_0^{4/9} B^{-0.842}$. Hence the ablation timescale $\tau_{abl} \simeq N_0/G_{tot} \propto N_0^{5/9} B^{0.842}$ goes up with the total injection amount and the magnetic field. If the injection amount is chosen to keep $N_0/V_{plasma}$ constant ($V_{plasma} \propto R^3$ is the plasma volume), then we have $\tau_{abl} \propto R^{5/3} B^{0.842}$. Thus even with the same thermal energy density and similar electron temperature, we would expect the larger (and those with stronger magnetic field) device plasma to be cooled down slower. Indeed, as is shown in Ref. [14], such effect could result in electron temperature collapse time on the order of 10ms in ITER. On top of such effect, the current contraction time in large devices is also inherently longer than that of small devices, due to the $L_\eta^2$ dependence in the resistive diffusion time where $L_\eta$ is the length scale of current redistribution. Meanwhile, the fragment flying time is only linearly depending on $L_\eta$. Hence even if the cooling rate is similar, the current redistribution time over the fragment flying time in large devices will still be larger than that of smaller devices, and relatively slower MHD response should be expected in larger devices as a result.

Both of the above two “size effects” would tend to result in slower current perturbation for larger devices. Such behaviour has raised some interest in exploring the possibility of temperature dilution without triggering the TQ promptly using D2 SPI [35]. Whether it is possible to utilize such effect to achieve sufficiently high radiation fraction while avoiding triggering the TQ prematurely is a interesting topic that will be further explored in future works.

### 3. The characteristic plasma response after impurity SPI

In this section, we demonstrate the characteristic dynamics after mixed impurity SPI into the “standard equilibrium” described in Section 2.2, using the parameters and numerical treatments described in Section 2.1 and Section 2.3. The MHD response caused by the injection induced cooling and its consequential density transport is described in Section 3.1, a discussion concerning the $n=0$ component of temperature evolution, as well as the transient 3D temperature structures and their impact is
presented in Section 3.2. Last, the existence of a stochastic transport barrier outside of
the $q = 1$ surface and its breaking down in the later phase of the TQ is discussed in
Section 3.3.

### 3.1. Characteristic MHD response and density transport

As has been explored in previous works [6, 13], the major MHD destabilization
mechanism after SPIs is the current perturbation caused by electron cooling, which
could be either axisymmetric, helical on major rational surfaces, or both. In the present
cases, it is found that due to the long current contraction time as a result of the “size
effects” discussed in Section 2.4, the $n = 0$ axisymmetric effect plays a minimal role in
the pre-TQ and early TQ phase.

In this subsection, we consider SPI with reference velocity 500m/s into the
“standard equilibrium” as an illustration of the characteristic MHD response and
the ensuing density transport. We calculate the NGS ablation rate using the local
plasma parameters, and the toroidal elongation parameter $\Delta\phi_{NG} = 0.5$. We call this
the “reference case”. The SPI simulations with the same equilibrium but different
injection velocities, ablation models and $\Delta\phi_{NG}$ show similar pre-TQ and early TQ
MHD dynamics, although with different timescale, different assimilation and radiation
strength. The electron and ion temperature cross-sections at the three times $t = 1.47$ms,
$t = 2.09$ms and $t = 2.45$ms are shown in Fig. 3, covering the pre-TQ and early TQ
phase. In agreement with previous Coronal Equilibrium (CE) impurity results [6], the
temperature profile shows a “scraping-off” behaviour as the fragments fly inward during
the pre-TQ phase, although the radiative cooling strength is stronger in our non-CE
case due to the larger fraction of low charge state neon [19], indicating a faster radiative
collapse onset. Then the core collapse is triggered as the fragments arrive on a major
rational surface in the core, in this case the $q = 1$ surface. The 1/1 kink goes unstable and
the ensuing flux displacement and “cold-bubble” [36] result in the temperature relaxation
within the $q = 1$ surface. The ion temperature is less relaxed than the electron one due
to the difference in their respective Braginskii parallel thermal conduction. Notably, the
temperature is not relaxed across the whole plasma, at least in the early TQ phase, as
can be seen from the evident ion temperature gradient shown in Fig. 3(f). This is due
to a transport barrier in the stochastic field which will be discussed later in this section.

The corresponding MHD spectrum for the first six toroidal harmonics is shown in
Fig. 4(a). For comparison, the MHD spectrum of a 200m/s case with the “non-local”
ablation scheme as introduced in Section 2.1 is also shown in Fig. 4(b). The solid lines are
the perturbed magnetic energies, and the dashed lines represent the perturbed kinetic
energies, while the colours represent different toroidal harmonics. The vertical chained
lines mark the time of strong MHD amplitude at the beginning of the core relaxation
and later on when the strong radiative collapse occurred. The first vertical chained line
in both cases corresponds to the vanguard fragments entering the $q = 1$ surface. Note
that there are two vertical lines in Fig. 4(b) simply because we run that case longer,
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The electron temperature cross-section at the same toroidal position as the SPI for the “reference case” at time (a) 1.47ms, (b) 2.09ms (the onset of the TQ) & (c) 2.45ms. The ion temperature cross-section at the same times with the electron ones are shown in (d)-(f). The SPI occurs at $t = 0$ms and the fragments are flying inward along the major radial direction from the LFS.</figcaption></figure>

and the temperature radiative collapse has the time to set in, resulting in the second MHD peak. In both cases the $n = 1$ modes are dominant. One notable feature here is that the 500m/s case shows weaker MHD amplitude compared with the 200m/s case when the core relaxation is triggered by fragments entering the $q = 1$ surface. This is consistent with the past numerical observation that slower injection results in faster radiative temperature collapse to 10eV range thus stronger current perturbation [14,15], although the whole pre-TQ phase of the slow injection case takes longer simply due to the longer fragment flying time.

As a side note, the plasma thermal energy is not depleted immediately as the core relaxation occurs due to the size-effect we discussed, but takes several milliseconds instead, as is shown in Fig. 5 for the 200m/s case with $\Delta\phi_{NG} = 0.3$. The comparison between the total thermal energy loss and the radiated energy, as well as the Ohmic heating is plotted. The TQ can be seen to occur gradually and in several stages, which is also reflected in the different MHD peaks in Fig. 4(b). We hereby define the beginning
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The perturbed energy spectrum after SPI for (a) the 500m/s local ablation case ("reference case") and (b) the 200m/s non-local ablation case. The solid lines represent the magnetic energies and the dashed lines are the kinetic energies, while different colours represent different toroidal harmonics. The vertical red chained line represent the time of MHD amplitude peaks.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The comparison between the total thermal energy loss (blue), the radiated energy (black), the Ohmic heating (red) and the combination of the two (magenta) for the 200m/s, nonlocal, $\Delta\phi_{NG} = 0.3$ case.</figcaption></figure>

of the core relaxation as the beginning of the TQ, but we remind the readers that it does not correspond to a fast depletion of all thermal energy. To discuss a bit more regarding the energy budget during the TQ, it can be seen that there was 6MJ of difference between the thermal energy loss and the radiated energy at $t \simeq 4.85$ms, after the core relaxation (not shown in the figure, the parallel kinetic energy and the ionization energy are both on the order of 1MJ at that time), resulting in a low radiated fraction. However, as we go deeper into the TQ, the radiated fraction rises as more energy is lost via radiation as can be seen by comparing the blue curve and the magenta curve, and the radiated fraction is about 80% by the end of the simulation. At this time, 25MJ out of the initial 36MJ thermal energy is lost.

The strong $n = 1$ MHD response is reflected in the density transport. The
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The electron density cross-section at the same toroidal position as the SPI for the “reference case” at times (a) 2.09ms, (b) 2.18ms (the onset of the TQ) & (c) 2.31ms.</figcaption></figure>

convective density transport at the time of the TQ onset is shown in Fig. 6. At the time of the TQ onset, significant density deposition exists on the $q = 1$ surface as marked by the density peak in Fig. 6(a). As the 1/1 kink develops, density cells can be seen “dragging” into the $q = 1$ surface as shown in Fig. 6(b) and (c). The density cells were transported over a distance of $O(1m)$ in 0.2ms, during which the fragments only travel 0.1m. Such MHD response and density transport behaviours are consistent with our previous observation with CE impurity models, where the $q = 1$ surface covers significant volume of the plasma [6].

To better demonstrate the core mixing effect caused by the kink motion, the comparison between the $n = 0$ component of the 500m/s non-local ablation 3D case electron density and that of an identical simulation but in 2D (assuming axisymmetry) is presented in Fig. 7. As would be expected, early in time due to the lack of MHD activity, there is not much difference between the two cases. As the fragments penetrate and MHD amplitude begins to rise, some inward density transport is seen in the 3D case as is shown by the green curve comparison between Fig. 7(a) & (b). As the 1/1 kink occurs, strong core density transport happens in the 3D case while the 2D case still shows a hollow electron density profile.

Overall, $n = 1$ modes are dominant in triggering the TQ onset, with significant 1/1 kink motion observed at the time of core relaxation. Such kink motion results in core density mixing and temperature relaxation. The electron temperature is more relaxed compared with the ion one due to the different parallel thermal conduction, but both experience a transport barrier on the $q = 1$ surface, which might be related with the incomplete destruction of good KAM surfaces [37–39] due to relatively small edge mode amplitude in the early TQ phase. The slower injection case with more significant radiative collapse [14, 15] shows stronger MHD amplitude at the time of core relaxation. This is partly due to the modes have more time to grow for the slower injection case, but more importantly the stronger radiative cooling results in stronger current perturbation
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Mid-plane cuts of the $n = 0$ axisymmetric electron density profiles for (a) the 500m/s, non-local ablation 3D case and (b) the 500m/s, non-local ablation 2D case. The different colours represent the profiles at different times.</figcaption></figure>

thus stronger drive for the instabilities.

## 3.2. The occurrence of radiative collapse and transient temperature structures

Since the current driven MHD amplitude is linked with the strength of radiative collapse, a closer look at the cooling dynamics and the transient temperature structure will offer more understanding on the MHD response. The electron cooling dynamics may also exhibit complicated interplay with the injection assimilation since the fragment ablation rate is also strongly depending on the ambient electron temperature [24].

As mentioned in the Section 1, it has already been found in axisymmetric simulations that, for fast enough fragment velocity and including the ablation reduction effect from the magnetic field as described by the $ζ_B$ multiplier in Eq. (14), ITER plasmas exhibit a delay of electron temperature radiative collapse after the passing of the vanguard fragments [14, 15]. Such a delay could be on the order of milliseconds, so that it is possible for the fragments to penetrate the plasma core before the edge temperature collapse.

On the other hand, 3D localization of the impurity density $n_{imp}$ and the electron density $n_e$ could have important impact on the radiative collapse, since the radiative power density is dependent on the local electron density according to Eq. (17), hence could exhibit strong peaking around the fragments in reality. It is therefore worthwhile to look into the details of the edge cooling in our 3D simulations and compare with the previous 1.5D and 2D results, and discuss their impact on the MHD response.

Mid-plane profiles of the $n = 0$ component of the ion and electron temperatures at different times for the reference case and a 200m/s non-local ablation case are shown in Fig. 8. The different colours represent profiles at different times, while the solid lines and dashed lines represent the electron and the ion temperature, respectively.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mid-plane cuts of the $n = 0$ axisymmetric electron (solid) and ion (dashed) temperature profiles for (a) the 500m/s, local ablation case and (b) the 200m/s, non-local ablation case. The different colours represent the profiles at different time.</figcaption></figure>

The most notable feature of the time evolution of both temperature profiles for the reference case is that the post-injection region remains at $O(100\text{eV})$ range as is shown in Fig. 8(a), consistent with the observed delay in the radiative collapse in previous lower dimension results [14, 15]. There is a transient increase of $T_e$ at the time of the red lines, corresponding to when the 1/1 kink begins to manifest. This is possibly a result of asymmetric flux tube compression [40] and high parallel electron thermal conduction which introduces the aforementioned compressional heating into the axisymmetric component. Accurately treating this process requires full MHD and also potentially inclusion of kinetic effects [40], which is beyond the scope of our paper, and we will leave it to future more detailed analysis with full MHD version of JOREK [18]. Another notable feature is the hollow temperature profile shortly after the onset of the TQ ($t = 2.12\text{ms}$), which is caused by the “cold-bubble” convecting the cold plasma into the core. As will be shown later in Section 3.3, most of the plasma is already stochastic at this time, and the electron temperature is apparently more relaxed than the ion one during the process, due to the higher parallel thermal conduction. Last, some transport barrier at $R = 5\text{m}$ and $R = 7.5\text{m}$ can be seen at $t = 2.12\text{ms}$ and $t = 2.45\text{ms}$ cases, roughly corresponding to the position of the $q = 1$ surface. This is the same transport barrier we mentioned in Section 3.1. As a comparison, the edge temperature plateau is less obvious in Fig. 8(b), and the radiative cooling in the outer region is more pronounced at the time of the core relaxation with almost no temperature plateau. This stronger radiation for the slow injection case can be corroborated by the total radiated energy in both cases as a function of the fragment penetration, which is simply estimated by multiplying the mean injection velocity with time, as is plotted in Fig. 9. Here, the radiated energy and neon assimilation for the above two cases and an axisymmetric 2D case are compared. While the neon assimilation is comparable for all three cases as
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The (a) radiated energy and (b) neon assimilation comparison between a 500m/s axisymmetric case (black), the 500m/s 3D case with local ablation scheme (blue) and the 200m/s 3D case with non-local ablation scheme (green).</figcaption></figure>

is shown in Fig. 9(b), the 200m/s case exhibits remarkably higher radiated energy for the same fragment penetration (thus cooling front position), consistent with the above temperature profile comparison as well as previous lower dimension results [14, 15].

The aforementioned delay in the edge radiative collapse, as well as the size-effect in the current contraction suggest the $n = 0$ current contraction is minimal [13]. Indeed it is the case as is shown in Fig. 10, where the toroidal current density $J_φ$ profile multiplied by the major radius $R$ is shown for both the 500m/s (reference case) and the 200m/s cases. Comparing Fig. 8(a) and Fig. 10(a), one may notice that despite the inward propagation of the cooling front from $t = 0$ms to $t = 1.47$ms, there is no axisymmetric current contraction, instead only some “jagged” current perturbations can be seen. These “jagged” current perturbations correspond to rational surfaces and indeed are consistent with our previous investigations in mildly cooling SPI cases [13]. Further, after the core relaxation occurred at $t = 2$ms, significant current flattening ensued within the $q = 1$ surface which is accompanied by a strong negative current sheet outside of the flattening region, as is represented by the magenta line. For the 200m/s case shown in Fig. 10(b), the edge $n = 0$ current is slightly more perturbed before $t = 4.81$ms due to the stronger radiative cooling, but still no significant $n = 0$ current contraction occurs due to the size effect in current redistribution. The current evolution dynamics looks similar to the 500m/s case, although the MHD amplitude is stronger for the 200m/s case. As the TQ proceeds, the radiative collapse and the current redistribution further develop, the edge current becomes more perturbed. Here one begins to see some “scraping-off” feature of the edge current density at $R = 4.5$m and $R = 8$m as anticipated from previous simplified model [41]. This corresponds to the second peak in the MHD amplitude shown in Fig. 4(b).

Despite the mild axisymmetric cooling, localized, transient electron temperature structures exist during the pre-TQ and early TQ phase. This is remarkable since the parallel electron thermal conduction is normally strong enough to relax such structures along the field lines. However, the strong local cooling and the relatively low thermal conduction in the high density, low temperature region around the fragments result in
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Mid-plane cuts of the $n = 0$ axisymmetric toroidal current $J_φ$ profile times the major radius $R$ for (a) the 500m/s, local ablation case and (b) the 200m/s, non-local ablation case. The different colours represent the profiles at different time.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The log$_{10}$ $T_e$ transient electron temperature structures at time $t = 1.26$ms for the 500m/s, non-local ablation case. (a) The Braginskii parallel thermal conduction case and (b) the constant parallel thermal conduction case which corresponds to that of the 500eV Braginskii parallel thermal conduction.</figcaption></figure>

local “electron temperature holes”. An example of such structure during the pre-TQ phase at $t = 1.26$ms for the 500m/s, non-local ablation case is demonstrated in Fig. 11 and Fig. 12.

In Fig. 11, a comparison of the log$_{10}$ $T_e$ between our Braginskii parallel thermal conduction case (Fig. 11(a)) and a constant thermal conduction case (Fig. 11(b)) at $t = 1.26$ms is shown, with all other parameters the same. The constant thermal conduction is equivalent to that of the 500eV Braginskii value. The most notable feature in the Braginskii case is a dark blue spot around $R \simeq 7.75m$ and $Z \simeq 0.65m$ in Fig. 11(a), which is absent in Fig. 11(b). Such a “electron temperature hole” corresponds to the
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The transient electron temperature and density structures at time $t = 1.26\,ms$ for the $500\,m/s$, non-local ablation case. (a) The electron temperature along several field lines started at various $R$, with the same $Z = 0.65\,m$, $\phi_{init} = \phi_{injector}$, as a function of the toroidal angle away from the injector; (b) The density along the field line for the same field lines.</figcaption></figure>

position of the fragment plume and is due to the reduced Braginskii thermal conduction in the cold region where a strong radiative heat sink exists.

A detailed look at the unrelaxed parallel electron temperature and density at the same time as Fig. 11 is shown in Fig. 12, where the electron temperature (Fig. 12(a)) and the electron density (Fig. 12(b)) are traced along field lines and shown as a function of the toroidal angle away from the initial position. The initial positions are chosen at different major radius, $Z = 0.65\,m$ and the same toroidal position as the fragment plume. Comparing the black and the red curves in Fig. 12(a), one can see the temperature hole develops as the electrons are rapidly cooling down and the Braginskii thermal conduction decreases. It further reaches very low temperature in the magenta curve, but as one moves further away from the vanguard fragments outward, the parallel temperature profile relaxes. This is due to the reduction of the radiative heat sink as the local ablation stops and the impurity density relaxes, hence the $T_e$ hole disappears after some time once the fragments pass.

Such strong parallel electron temperature gradient naturally raises the question of whether the collisional Braginskii thermal conduction is still valid and whether the electron distribution within the temperature hole is close to Maxwellian. To resolve this question, the inverse of the hot electron collisional cross-section could be used to estimate the required line integrated density for its thermalization. We consider the electron temperature far from the temperature hole to be $T_e = 500\,eV$, corresponding to the situation of the red curve in Fig. 12(a). Any cases colder than the red curve case should be more collisional so that the Braginskii description should be more valid there. The required line integrated density for a collision to occur is

$$v_{te}(T_e)n_e/\nu_{ee}(n_e, T_e) = \langle\sigma\rangle^{-1} \approx 2 \times 10^{21}/m^2. \qquad (20)$$

Here $v_{te}(T_e)$ is the ambient electron thermal velocity. This requirement should be
compared against integration along field line across the half-width of the temperature hole. Here, we roughly estimate it by taking an integration along the red curve of Fig. 12(b) from $\delta\phi = 0$ to $\delta\phi = 1$, approximating the major radius as $R \approx 6m$. The resulting line integrated density is $\int n_e dl \approx 2.4 \times 10^{21}/m^2$, comparable with the $(\sigma^{-1})$ obtained above. Thus the electron heat transport into the temperature hole is barely collisional even when they are still relatively warm, so that we still consider our use of the Braginskii thermal conduction a good enough description of the electron parallel heat transport. On the other hand, past studies have pointed out that the hot electron thermalization with the bulk, cold electrons could be approximated as the rapid decay of the hot species density while its mean energy remains unchanged or even increases [42–44]. Hence we expect the electron distribution within the electron temperature hole to be mostly Maxwellian, with a small long-tail population. These long-tail electrons could penetrate the plasmoid around the fragments and significantly contribute to the ablation process [25,26]. To take this effect into account, we introduced the “non-local” ablation scheme in Section 2.1 so that the flux averaged electron parameters are used in calculating the ablation rate, instead of the local parameters.

It is also noteworthy that the electron temperature hole shown in Fig. 11 and Fig. 12 will cause enhanced current expulsion from the field lines they are on, as they enhance the resistivity along the field line. This causes the Braginskii cases to experience relatively stronger MHD amplitude compared with the constant thermal conduction cases as will be seen in Section 4.4.

### 3.3. The stochastic transport barrier on the $q=1$ surface

The last important aspect of this section is the transport barrier outside of the $q=1$ surface observed in Fig. 3 and Fig. 8. Indeed, we will see that this is a reflection of the corresponding transport barrier in the stochastic field lines. To demonstrate this, we initialize field line tracers evenly positioned from the magnetic axis to the plasma boundary both for our 500m/s “reference” case studied in this section, and also for the SPI case with 200m/s injection velocity and non-local ablation scheme. The MHD spectrum of both cases have been shown in Fig. 4. We trace the field lines along both directions for 2000 toroidal turns and plot the resulting Poincaré plot for the above two cases in Fig. 13, with the point colour representing different field line tracers in Fig. 13(a), (b) &(c), and different remaining parallel length until the field line hits the wall in Fig. 13(d), (e) &(f). If the field line remains well confined after 2000 turns, it is marked red in Fig. 13(d), (e) &(f). Here Fig. 13(a) is the 500m/s, local ablation scheme case with $\Delta\phi_{NG} = 0.5$ at $t = 2.39$ms, some time after the core relaxation. It can be seen that although the field lines are already stochastic, the field lines within the $q=1$ surface do not easily penetrate outside as is indicated by the concentration of purple colours within. The connection length plotted in Fig. 13(d) shows consistent feature with (a), as the field lines within the $q=1$ surface are all well confined after 2000 turns. This behaviour is reflected in the transport barrier shown in Fig. 3 and Fig. 8. This
<figure><img src="image_14.png" /><figcaption>**Figure 13.** The Poincaré plot for (a) & (d) the 500m/s, local ablation scheme case with $\Delta\phi_{NG} = 0.5$ at time $t = 2.39\,\text{ms}$; (b) & (e) the 200m/s, non-local ablation scheme case with $\Delta\phi_{NG} = 0.3$ at time $t = 4.85\,\text{ms}$; and (c) & (f) the 200m/s, non-local ablation scheme case with $\Delta\phi_{NG} = 0.3$ at time $t = 6.36\,\text{ms}$. The colours in (a), (b), (c) represent different field line tracers. The colours in (d), (e), (f) represent different remaining field line length until the field line hits the wall, while the well-confined field lines after 2000 toroidal turns are marked as red.</figcaption></figure>

is likely caused by the stiff dominant 1/1 kink mode structure, which quickly vanishes outside of the $q = 1$ surface. Meanwhile, the edge mode amplitude, such as the 2/1 or 3/1 mode amplitude, is weak due to the mild edge cooling discussed above. Both factors combine to result in relatively small magnetic perturbation outside of the $q = 1$ surface, so that some good KAM surfaces [37–39] might survive, limiting the outward penetration of the stochastic field lines.

Such transport barrier is a common feature for our “standard equilibrium” which has a large $q < 1$ region. In Fig. 13(b), the same plot is shown for the 200m/s, non-local
ablation scheme case with $\Delta\phi_{NG} = 0.3$ at $t = 4.85\,\text{ms}$, shortly after the core relaxation. The edge field lines shown in Fig. 13(b) & (e) are more stochastic than those of Fig. 13(a) & (d), as is suggested by their better mixing and shorter connection length. This is due to the stronger MHD amplitude as is shown in Fig. 4(b). Despite this, a similar stochastic transport barrier is still present around the $q = 1$ surface, as is demonstrated by the high concentration of blue and purple Poincaré points within said surface in Fig. 13(b) and the persistent, although reduced, well-confined region in Fig. 13(e). As the TQ proceeds, however, further growth of the MHD amplitude caused by the continued edge cooling finally breaks the transport barrier as is shown in Fig. 13(c) & (f), representing the $200\,\text{m/s}$, non-local ablation scheme case at $t = 6.36\,\text{ms}$. There, good mixing between different tracers can be seen across the plasma volume, and no field line remains well-confined. No particular concentration of Poincaré points is observed, apart from some small remnant flux tubes, and the connection length to the wall is on the order of $10^4\,\text{m}$. Despite not running long enough in the $500\,\text{m/s}$ case, we anticipate to see a similar breaking-down of the transport barrier after the periphery plasma is sufficiently cooled down, although this may take on the order of $10\,\text{ms}$ [14].

To further demonstrate the existence and breaking-down of this transport barrier, we initialize 300 field line tracers evenly positioned along the toroidal direction, with initial poloidal location $R_{init} = 6.4\,\text{m}$ or $R_{init} = 7.5\,\text{m}$ and $Z_{init} = 0.6\,\text{m}$ for the $200\,\text{m/s}$ case discussed above. The resultant Poincaré plot is shown in $(\theta, \psi_n^{1/2})$ space in Fig. 14, where $\psi_n$ is the normalized flux and $\theta$ is the polar angle. Different colours here represent different line tracers. Fig. 14(a) and (c) show the plot for tracers originating from the core ($R_{init} = 6.4\,\text{m}$) at time $t = 4.85\,\text{ms}$ and $t = 6.05\,\text{ms}$ respectively, corresponding to the first and second MHD peaks in Fig. 4(b), while Fig. 14(b) and (d) correspond to the plots for tracers launched outside of the $q = 1$ surface at $R_{init} = 7.5\,\text{m}$. All tracers are followed for 500 toroidal turns. A clear separation at the black vertical chained line (the $q = 1$ surface) can be seen for Fig. 14(a) and (b), marking the transport barrier at $t = 4.85\,\text{ms}$. On the other hand, Fig. 14(c) and (d) show good mixing across the whole plasma volume, signifying the breaking-down of the barrier as the edge modes develop, which is also shown in Fig. 4(b) by the growing $n = 1$ and $n = 2$ modes. The change of $q = 1$ surface position is the combined effect of slight central current profile changes and the relatively flat $q$ profile in the plasma core.

The level of stochasticity can also be assessed by looking at the mean position and position variance of all the field line tracers as a function of toroidal turns as is shown in Fig. 15. The mean position and its variance for tracers launched from the core and the edge at $t = 4.85\,\text{ms}$ are shown in Fig. 15(a) and (b), while those at time $t = 6.05\,\text{ms}$ are shown in Fig. 15(c) and (d) respectively. We are only counting the field lines remaining in the plasma after 500 toroidal turns, and different colours represent different tracers. The transport barrier at $t = 4.85\,\text{ms}$ is evident in Fig. 15(a) as the blue and red curves are converging to different mean positions. The core stochasticity is stronger since the blue curve reaches its steady result within 100 toroidal turns in both Fig. 15(a) and (b), while it takes hundreds of toroidal turns for the red curve to spread out as indicated by
<figure><img src="image_15.png" /><figcaption>**Figure 15.** The mean position and its variance in the square root of normalized flux $\psi_n^{1/2}$ for the 200m/s case. The field line tracers are initially positioned equally along the toroidal direction, launched at $Z_{init} = 0.6m$ at time (a) & (b) $t = 4.85ms$ and (c) & (d) $t = 6.05ms$. Different colours represent different line tracers.</figcaption></figure>

## 4. Impact from different ablation model, toroidal localization, background impurities and parallel thermal conduction

Since the electron cooling dynamics are closely linked to the current perturbation which in turn determines the MHD amplitude and the consequent field line stochasticity, it is of interest to see how the various physical assumptions and numerical treatments affect the radiative power and thus the edge cooling. In Section 4.1, we shall investigate the impact from two ablation model assumptions, the first is the use of the flux-surface-averaged ambient plasma parameters instead of the local ones in calculating the ablation rate, the second is to include the reduction in the ablation rate due to the large magnetic field introduced in Section 2.1. Further, in Section 4.2, the assimilation and radiation for two different toroidal deposition widths are compared. The background impurity contribution in the radiative power is compared with that of the total radiative power in Section 4.3. Last, the use of Braginskii parallel thermal conduction is compared against that of a constant thermal conduction in Section 4.4.
4.1. *Impact from the non-local ablation and the magnetic effect*

We have discussed the transient “electron temperature hole” structure around the
fragments and its implication in Section 3.2, and introduced the non-local ablation
scheme in Section 2.1. Here, we compare the assimilation and radiation using the local
ablation scheme and the non-local one. The comparison of radiation power and the
total neon assimilation between the 3D local ablation scheme case, the 3D non-local
ablation scheme case and an axisymmetric case are shown in Fig. 16. All cases share the
same velocity spread and the same 500m/s reference injection velocity, and the two 3D
cases share the same density toroidal deposition width Δφ<sub>NG</sub> = 0.5. It can be seen from
Fig. 16(a) that the 3D local ablation case experiences similar radiation power compared
with the 2D case, while the 3D non-local ablation case sees a significant increase in the
radiation power. The reason for such behaviour is the following. On one hand, the
radiation power for each neon atom or ion is proportional to the local electron density,
such that 3D cases tend to experience more radiation due to their un-relaxed density
peak. On the other hand, the 3D local ablation case experience smaller ablation rate
due to the “electron temperature hole”, which is shown in Fig. 16(b), so that there are
less neon ions or atoms to contribute to the total radiation power. The two effects
approximately cancel each other in the early phase before t = 2.5ms, resulting in the
similar radiation power between the blue and the black lines. Meanwhile, the 3D non-
local ablation scheme experiences similar assimilation with the 2D case as is shown in
Fig. 16(b), so that only the density localization effect remains in the total radiation
power, resulting in the stronger radiation for the red curve shown in Fig. 16(a). Apart
from such localization effect, there could also be some contribution from the fact that
the 3D case is accompanied by field line stochasticity thus an outgoing heat flux into
the radiative region, thus increasing the radiation power. Indeed, it can be seen that
the radiation spikes at t = 1.75ms and t = 2.1ms in Fig. 16 roughly corresponds to the
MHD activities at the same time in Fig. 4. We do not see significant difference in the
ionization energy in general due to the robustness in the neon assimilation shown in
Fig. 16(b).

Comparing Fig. 16(a) and (b), one may find it counter-intuitive that the impact
of the non-local ablation scheme is more pronounced in the radiation power compared
with the assimilation. Such behaviour can be understood as follows. The ablation
process is self-regulated in the sense that stronger ablation will lower the ambient plasma
temperature more quickly and reduce the subsequent ablation. Hence the enhancement
caused by the non-local calculation did not have such a strong impact on the assimilation
itself, and indeed the non-local assimilation converges with that from the 2D case which
is as expected. On the other hand, the radiative power is not self-limiting before the
temperature radiative collapse occurs, hence in the ITER plasma we considered here the
impact on the radiative power is more apparent compared with that on the assimilation.

Although we are mainly concerned with the pre-TQ and early TQ phase in this
paper, it might be beneficial to also take a look at the ratio between the radiated
<figure><img src="image_16.png" /><figcaption>**Figure 16.** The (a) radiative power and (b) neon assimilation comparison between an axisymmetric case (black), a 3D case with local ablation scheme (blue) and a 3D case with non-local ablation scheme (red).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 17.** The comparison between the total thermal energy loss (chained) and the radiated energy (solid). The different colours represent the same cases as is shown in Fig. 16.</figcaption></figure>

energy and the total thermal energy loss, as is shown in Fig. 17. There, different colours represent the same cases as is shown in Fig. 16, and the solid curves are the radiated energy and the chained ones are the thermal energy loss. It is shown that while the total thermal energy loss is comparable between the local and the non-local cases, the increased radiation power for the non-local case means better radiated fraction. It should be borne in mind that at this early phase this fraction could be small compared with that by the end of the TQ, since most significant heat flux is expected when the relaxation first occurs. Indeed, as is shown previously in Fig. 5, as time goes on this radiated power fraction increases.

We now move on to investigate the effect of the magnetic constraining term $\zeta_B$ in Eq. 14. Two cases with and without such magnetic effect are compared in Fig. 18. Both cases have the same injection velocity 200m/s, both use the non-local ablation scheme and both have the same $\Delta\phi_{NG} = 0.3$. It can be seen in Fig. 18(a) that the case without the magnetic effect experiences stronger radiative power compared with
<figure><img src="image_18.png" /><figcaption>**Figure 18.** The (a) radiative power and (b) neon assimilation comparison between the 200m/s, non-local ablation scheme case with magnetic effect (black) and without (blue). Both cases have the same $\Delta\phi_{NG} = 0.3$.</figcaption></figure>

the case with the magnetic effect, while the former case also shows some enhancement in the assimilation compared with the latter case. The enhancement in the assimilation is, however, much weaker than one would have expected from the $\zeta_B \sim (\frac{2}{\max(B[T],2)})^{0.843}$ dependence. Such behaviour is consistent with our previous discussion concerning Fig. 16, that enhancement in the ablation rate would lead to only limited improvement in the assimilation due to its self-regulated nature, while the impact on the radiation could be more pronounced.

## 4.2. Impact from toroidal localization in density deposition

Since the radiative power for each impurity atom or ion is proportional to its local electron density, concerns naturally arise regarding the accuracy of radiation power produced by the mainstream 3D fluid codes since the density source deposition is usually artificially elongated along the toroidal direction due to numerical limitations in the toroidal resolution. To demonstrate the impact of such toroidal localization, we compare a 500m/s, 3D non-local ablation case with $\Delta\phi_{NG} = 0.3$ against the three cases discussed in Fig. 16 where $\Delta\phi_{NG} = 0.5$. The result is shown in Fig. 19.

As expected, the enhanced toroidal localization indeed causes stronger radiation power as can be seen by comparing the magenta and red curves in Fig. 19(a). The earlier peak in the magenta curve around $t = 1.3$ ms is due to an early MHD burst which causes some outward heat flux into the outer radiative region. The enhanced radiation in the magenta curve is purely due to enhanced density localization, as the total neon assimilation is similar between the red and the magenta curves as can be seen in Fig. 19(b). Meanwhile, the radiation power increase is less than proportional to the change in the $\Delta\phi_{NG}$. This is likely the result of competition between the local source and the parallel relaxation, which causes the density peak to be less localized than the shape of the density source. This could also be caused by the merging of the plasmoid from each fragment to form an envelop of a “plasmoid conglomerate” which smooth over the density deposition.
<figure><img src="image_19.png" /><figcaption>**Figure 19.** The (a) radiative power and (b) neon assimilation comparison for 500m/s cases. The black, blue and red curves represent the 2D, 3D local ablation and 3D non-local ablation cases plotted in Fig. 16. All three cases are with $\Delta\phi_{NG} = 0.5$. The new magenta curves represent a 3D non-local ablation case with $\Delta\phi_{NG} = 0.3$.</figcaption></figure>

Overall, we found that increased toroidal localization indeed results in enhanced radiation power, thus present and future fluid numerical radiation results with an artificial toroidal elongation should be considered as a lower estimation of the actual radiative power. On the other hand, the radiation increase is not linearly dependent on the toroidal localization of the density source, but shows a weaker dependence. Exploration on how to account for such localization effect given the limited toroidal resolution is a challenging task for future works.

### 4.3. Impact from the background impurities

Apart from the massive amount of impurities injected into the plasmas, the already existing background impurities may also contribute to the temperature radiative collapse, and it is of interest to quantify the relative importance of such background impurity effect. For ITER, we take a simplified model of uniform background impurity distribution. Two background impurities are considered, with beryllium density $n_{Be} = 3.5 \times 10^{17}/m^3$ and tungsten density $n_{W} = 3.5 \times 10^{11}/m^3$, according to the expectation of such L-mode discharge [29]. The background impurities are modelled using the Coronal Equilibrium assumption. Comparisons between the total impurity radiation power and that from the background species alone are shown in Fig. 20 for the 500m/s case and the 200m/s case respectively.

It is found that after the fragments enter the separatrix, the total radiation power rises sharply in both cases, but the contribution from the background species only rises very slowly, and there are approximately three orders of magnitude difference between the total radiation power and the background impurity contribution. Not shown in the figure, the majority of the background impurity radiation is coming from the Beryllium due to its relatively high density, and tungsten only accounts for 0.1 percent of the background impurity radiation. Although we consider a uniform background impurity distribution here, a non-uniform impurity distribution keeping
<figure><img src="image_20.png" /><figcaption>**Figure 20.** The comparison between the total radiative power from all impurity sources (black) and that from the background impurities alone (red) for (a) the 500m/s case and (b) the 200m/s case. All cases use the same $\Delta\phi_{NG} = 0.3$ and non-local ablation scheme.</figcaption></figure>

the same total impurity quantity would not impact our conclusion due to the large difference in population compared with the injected neon. This may change in the H mode case where tungsten density is on the order of $10^{15}/m^3$ for ITER. But it is still unlikely to be comparable with the radiation from the injected impurities. Thus we conclude that the background impurity radiation is a negligible effect in the occurrence of the electron temperature radiative collapse unless some events introduce drastically enhanced background impurity accumulation in the plasma.

### 4.4. Impact from the parallel thermal conduction

The last part of interest is how the parallel thermal conduction model affects the radiative cooling and assimilation. As has been shown in Fig. 11, constant parallel thermal conduction prevents the formation of significant electron temperature holes. This means the perturbation to the parallel current is less significant in the constant thermal conduction case, thus the MHD response is weaker. To demonstrate this, we show the Poincaré plot of a 500m/s, non-local ablation, $\Delta\phi_{NG} = 0.3$ case with constant thermal conduction corresponding to that of the Braginskii value at $T_e = 500$eV in Fig. 21. The time is $t = 2.39$ms, after the vanguard fragments have entered the $q=1$ surface and destabilized the 1/1 kink. The field line colours here represent different initial positions of field line tracers. It can be seen that although the flux surfaces within the $q=1$ surface are destroyed by the 1/1 kink, relatively good flux surfaces remain in the outside region, as opposed to the situation in Fig. 13(a) where, despite the $q=1$ transport barrier, the outer region is already stochastic.

The impact of such weakened MHD response can be seen in Fig. 22, which is the same as Fig. 19 but with an additional cyan curve that represents our constant thermal conduction case. Apart from the similar assimilation shown in Fig. 22(b), the constant thermal conduction case shows a delayed rise of the radiation power as can be seen in Fig. 22(a). Before $t = 2$ms, the constant thermal conduction case shows weaker radiation power due to the lack of MHD induced outward heat flux from the core, while
<figure><img src="image_21.png" /></figure>

**Figure 21.** The Poincaré plot at $t = 2.39\text{ms}$ for the constant thermal conduction case, after the fragments enter the $q = 1$ surface. Different colours represent difference field line tracers.

<figure><img src="image_22.png" /></figure>

**Figure 22.** The (a) radiative power and (b) neon assimilation comparison for $500\text{m/s}$ cases. The black, blue, red and magenta curves represent the same cases plotted in Fig. 19. The new cyan curves represents a 3D non-local ablation case with $\Delta\phi_{NG} = 0.3$ and constant thermal conduction corresponding to that of Braginskii at $500\text{eV}$.

the other 3D non-local cases already show burst like radiation peaks. Only after the fragments entered the $q = 1$ surface and triggered the 1/1 kink does one see the radiation peak caused by the core relaxation which is on similar level with the Braginskii case.

## 5. Conclusion and discussion

In conclusion, emerging studies have point to the possibility of fragments penetrating faster than the inward propagation of the temperature radiative collapse in large machines such as ITER. Since violent MHD instabilities which result in strong stochastic
transport during the TQ phase are usually triggered by the radiative collapse on major rational surfaces, the above behaviour point to the possibility of dilution cooling of the bulk plasma before these violent MHD response occurs, thus helping in reducing the hot-tail runaway electron seed as well as helping to reduce the TQ thermal load. Some important 3D features such as the localization effect and the MHD induced transport are overlooked in previous lower dimension works, however, motivating us to look into the the neon assimilation, the radiative cooling and the accompanying MHD response in the pre-TQ and early TQ phase of SPI in ITER plasmas with 3D simulation. Our result confirms that the fragment can indeed reach the core before strong MHD responses and stochastic transport are triggered by the radiative cold front, especially for the 500m/s fast injection case. The characteristic post-injection dynamics are found to be quite similar despite differences in the injection parameters, physical assumptions and numerical treatment, although the exact cooling strength and the corresponding MHD amplitude are different.

The dominant mode is $n = 1$ as expected for single SPI cases. Due to the presence of a large $q < 1$ region, the core temperature relaxation and the dominantly convective density transport show features distinctive of a 1/1 kink which is triggered by the vanguard fragments entering the corresponding surface, in agreement with past simulations with Coronal Equilibrium assumptions [6, 13]. The MHD amplitude during such relaxation generally depends on the radiative cooling strength, with stronger cooling causing stronger current perturbation thus stronger MHD amplitude. Due to the size effect in large devices, the axisymmetric cooling (thus the axisymmetric current contraction) is mild in most cases investigated, with faster injections showing less cooling and MHD response compared to the slower ones with the same fragment penetration. Such behaviour is in agreement with previous 1.5D and 2D results [14]. However, the electron temperature profile also shows significant transient, 3D “electron temperature hole” structures. Such structures have two effects, one is resulting in stronger current perturbation due to the increased resistivity in the parallel current path, thus stronger MHD response; the other is motivating the non-local consideration of the fragment ablation calculation. Another notable feature is the transport barrier in stochastic field lines outside of the $q = 1$ surface, which originates from the combination of the stiff 1/1 mode structure and the weak edge mode amplitude as a result of the aforementioned mild cooling. As time goes on and the temperature radiative collapse occurs, the edge MHD amplitude ultimately grows enough to destroy such barrier and cause good mixing of the stochastic field lines. This occurs towards the late phase of the TQ as is shown by combining Fig. 4(b), Fig. 13 and Fig. 5, at the time of the transport barrier breaking-down after $t = 6$ms by the second MHD peak, 2/3 of the initial thermal energy is already radiated away.

Four different factors are considered for their impact on the radiative cooling strength in the pre-TQ and early TQ phase. The ablation model, the toroidal density deposition width, the background impurity contribution and the parallel thermal conduction model. For the ablation model, it is found that the use of the non-
local ablation scheme and the exclusion of the magnetic effect result in only a mild increase of the assimilation due to the self-regulated nature of the ablation. Meanwhile, the radiation power shows stronger increase since it is not self-regulated before the occurrence of the temperature radiative collapse. The more localized toroidal density deposition is found to indeed result in enhanced radiation due to its dependence on the electron density. Such enhancement is not proportional to the density source localization, however, probably due to the density peaking is mitigated by the parallel relaxation or the merging of the plasmoid clouds. The background impurity radiation is found to be negligible compared with that from the injected impurities due to the huge difference in density. The use of constant parallel thermal conduction is found to mitigate the occurrence of “electron temperature hole” seen in the Braginskii case, leading to a weaker MHD response and weaker stochasticity.

Based on above numerical results, it appears that the occurrence of the temperature collapse and strong current perturbation is facilitated by any enhancement in the ablation rate, be it the exclusion of the magnetic effect, the non-local scheme or simply increasing the fragment surface-to-volume ratio. The other way of facilitating the radiative collapse is the consideration of realistic toroidal elongation of the density source, which is out of reach for major MHD codes due to limits regarding the toroidal resolution at present. The radiation prediction of any simulation with elongated toroidal density deposition should be considered as a lower estimation compared with reality, although the radiation enhancement is less than linearly proportional to the density localization as discussed in Section 4.2. How best to account for the radiation power in reality with currently limited toroidal resolution is a concern for future works.

There are other factors which could contribute to the occurrence of the radiative collapse but are not considered in our current study. One is the presence of gaseous materials or small droplets after the pellet shattering, which could travel faster than the bulk fragment plume and result in stronger edge cooling which facilitate the edge radiative collapse. The other is the dependence of the characteristic fragment size on the injection velocity, which tends to produce finer fragments for faster injections. This may result in larger surface-to-volume ratio for the faster injection cases thus countering the mild cooling behaviour we observed here for the fast injection cases. One could usually modify the shattering angle to mitigate such enhancement, however. Another factor is the $q$ profile of the target plasma, in the absence of the large $q = 1$ mode one might anticipate that the 1/1 will not play any important role, but 2/1 or 3/2 would dominate instead. With their different mode structures compare with that of the 1/1, one might expect to see different transport behaviour of the stochastic field lines as discussed in this study. For example, one would anticipate relatively weaker core density transport at the time of the core relaxation due to the absence of the strong kink motion, but the stochastic field lines may have an easier time to penetrate outside of the resonant surface compared with the 1/1 dominant cases since the 2/1 and 3/2 mode structure are not as stiff. Last, an injection velocity on the order of 1000m/s might be achievable with future electro-magnetic injectors, so that the fragments, or
cluster of small pellets, may penetrate into the high field side before the periphery plasma experiences radiative temperature collapse, results in further cooling down of the periphery plasma and stronger MHD amplitude. This may counteract the milder cooling effect of faster injections. The inclusion of those further effects would be pursued in future works.

Extrapolating to the ITER H-mode SPI scenario, we expect the higher thermal energy density in H mode to cause enhanced assimilation and stronger un-relaxed density peak around the fragments, thus stronger local radiation peak due its dependence on the electron density. In the case of very small neon mixture ratio, if the local radiation is not enough to deplete the strong thermal reservoir quickly, then we would expect to see strong local pressure peak which may drive plasmoid drift outward along the major radial direction [45], which may decrease the neon assimilation within the separatrix. This effect should tend to vanish in the large neon mixture-ratio limit due to the reduction of the local pressure peak by radiation [46]. Further, we also expect the H mode pedestal to be strongly destabilized by the vanguard fragments [47], this would likely to cause premature thermal energy loss before deep fragment penetration. On the other hand, whether we will really perform SPIs into a full H-model in reality is uncertain, as significant thermal energy loss would likely occur already in the pre-TQ phase of disruptions, and many disruptive discharges would see the transition back to L-mode before the SPI system is triggered [48]. Hence a more realistic scenario would be injection into a post-H-mode case where the pedestal is already collapsed. Exploring the H-mode SPI scenario is an immediate future work we are going to pursue.

Acknowledgments

The authors thank A. Matsuyama, P. Aleynikov, C. C. Kim and M. Kong for fruitful discussion. ITER is the Nuclear Facility INB no. 174. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization. This publication is provided for scientific purposes only. Its contents should not be considered as commitments from the ITER Organization as a nuclear operator in the frame of the licensing process. Part of this work is supported by the National MCF Energy R&D Program of China under Grant No. 2019YFE03010001 and the National Natural Science Foundation of China under Grant No. 11905004. This work has been co-funded by the ITER Organization under the implementing agreement IO/IA/19/4300002053. Part of this work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement 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. This work is carried out on Tianhe-3F operated by NSCC-TJ.

References
[1] N. Commaux, D. Shiraki, L.R. Baylor, et al., "First demonstration of rapid shutdown using neon shattered pellet injection for thermal quench mitigation on DIII-D", Nucl. Fusion 56 046007 2016;

[2] D. Shiraki, N. Commaux, L. R. Baylor et al., "Thermal quench mitigation and current quench control by injection of mixed species shattered pellets in DIII-D", Phys. Plasmas 23 062516 (2016);

[3] Y. Li, Z.Y. Chen, W. Yan et al., "Comparison of disruption mitigation from shattered pellet injection with massive gas injection on J-TEXT", Nucl. Fusion 61 126025 (2020).

[4] S. Jachmich, U. Kruezi, M. Lehnen et al., "Shattered pellet injection experiments at JET in support of the ITER disruption mitigation system design", Nucl. Fusion 62 026012 (2022).

[5] C. C. Kim, Y. Liu, P. B. Parks et al., "Shattered pellet injection simulations with NIMROD" Phys. Plasmas, 26, 042510 (2019).

[6] D. Hu, E. Nardon, M. Hoelzl et al., "Radiation asymmetry and MHD destabilization during the thermal quench after impurity Shattered Pellet Injection" Nucl. Fusion, 61, 026015 (2021).

[7] J. McClenaghan, B.C Lyons, C.C. Kim et al., "MHD modeling of shattered pellet injection in JET" Nucl. Fusion, submitted (2023).

[8] E. M. Hollmann, T. C. Jernigan, E. J. Strait, et al., "Observation of q-profile dependence in noble gas injection radiative shutdown times in DIII-D", Phys. Plasmas 14, 012502, 2007;

[9] A.J. Thornton, K.J. Gibson, I.T. Chapman, et al., "Plasma profile evolution during disruption mitigation via massive gas injection on MAST", Nucl. Fusion 52, 063018, 2012;

[10] D. Hu and L. E. Zakharov, "Quasilinear perturbed equilibria of resistive unstable current carrying plasma", J. Plasma Phys. 81 515810602 (2015);

[11] R. B. White, D. A. Gates and D. P. Brennan, "Thermal island destabilization and Greenwald density limit", Phys. Plasma 22 022514 (2015);

[12] Q. Teng, N. Ferraro, D.A. Gates et al., "Nonlinear simulations of thermo-resistive tearing mode formalism of the density limit", Nucl. Fusion 58 106024 (2018);

[13] D. Hu, E. Nardon, M. Lehnen et al., "3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection" Nucl. Fusion, 58, 126025 (2018).

[14] A. Matsuyama, D. Hu, M. Lehnen et al., "Transport simulations of pre-thermal quench shattered pellet injection in ITER: code verification and assessment of key trends" Plasma Phys. Control. Fusion 64, 105018 (2022).

[15] A. Matsuyama, "Cryogenic Pellet Ablation Physics and Integrated Modelling of Shattered Pellet Injection", 2nd IAEA Technical Meeting on Plasma Disruptions and Their Mitigation, 2022.

[16] P. Aleynikov, B. N. Breizman, Per Helander, et al., "Plasma ion heating by cryogenic pellet injection", J. Plasma Phys., vol. 85, 905850105 (2019).

[17] A. M Arnold, P. Aleynikov and P. Helander, "Self-similar expansion of a plasmoid supplied by pellet ablation", Plasma Phys. Control. Fusion 63 095008 (2021).

[18] M. Hoelzl, G.T.A. Huijsmans, S. Pamela et al., "The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas" Nucl. Fusion 61 065001 (2021).

[19] D Hu, G T A Huijsmans, E Nardon et al., "Collisional-radiative non-equilibrium impurity treatment for JOREK simulations" Plasma Phys. Control. Fusion 63 125003 (2021).

[20] H P Summers, W J Dickson, M G O'Mullane et al., "Ionization state, excited populations and emission of impurities in dynamic finite density plasmas: I. The generalized collisional–radiative model for light elements" Plasma Phys. Control. Fusion 48 263-293 (2006).

[21] H.P. Summers and R.W.P. McWhirter, "Radiative power loss from laboratory and astrophysical plasmas. I. Power loss from plasmas in steady-state ionisation balance", J. Phys. B: At. Mol. Phys. 12 2387, (1979);

[22] L. Spitzer, Physics of Fully Ionized Gases (Interscience, New York, 1956).

[23] S.P. Hirshman, "Neoclassical current in a toroidally-confined multispecies plasma", Phys. Fluids 21, 1295 (1978).
[24] J. Zhang and P. B. Parks, "Analytical formula for pellet fuel source density in toroidal plasma configurations based on an areal deposition model", Nucl. Fusion **60** 066027 (2020);

[25] P.B. Parks, R.J. Turnbull and C.A. Foster, "A model for the ablation rate of a solid hydrogen pellet in a plasma", Nucl. Fusion **17** 539 (1977).

[26] B Pégourié, V Waller, R J Dumont et al., "Modelling of pellet ablation in additionally heated plasmas" Plasma Phys. Control. Fusion **47** 17 (2005).

[27] S. I. Braginskii, "Transport Processes in a Plasma", Rev. Plasma Phys. **1**, 205, (1965);

[28] H.P. Summers, "Atomic data and analysis structure users manual, JET-IR 06, Abingdon: JET Joint undertaking", 1994. http://open.adas.ac.uk/

[29] S.H. Kim, T.A. Casper and J.A. Snipes, "Investigation of key parameters for the development of reliable ITER baseline operation scenarios using CORSICA", Nucl. Fusion **58** 056013 (2018);

[30] J P Graves, D Zullino, D Brunetti et al., "Reduced models for parallel magnetic field fluctuations and their impact on pressure gradient driven MHD instabilities in axisymmetric toroidal plasmas", Plasma Phys. Control. Fusion **61** 104003 (2019);

[31] F. Villone, Y.Q. Liu, G. Rubinacci et al., "Effects of thick blanket modules on the resistive wall modes stability in ITER", Nucl. Fusion **50** 125011 (2010);

[32] T. E. Gebhart, L. R. Baylor, and S. J. Meitner, "Analysis of the Shattered Pellet Injection Fragment Plumes Generated by Machine Specific Shatter Tube Designs", Fusion Sci. Technol., **77:1**, 33-41 (2021);

[33] P. B. Parks, "Modeling dynamics Fracture of Cryogenic Pellets", GA report GA-A28352, (2016).

[34] T.E. Gebhart, L.R. Baylor & S.J. Meitner, "Shatter Thresholds and Fragment Size Distributions of Deuterium-Neon Mixture Cryogenic Pellets for Tokamak Thermal Mitigation", Fusion Sci. Technol., 76:7 831-835 (2020).

[35] E. Nardon, D. Hu, M. Hoelzl et al., "Fast plasma dilution in ITER with pure deuterium shattered pellet injection", Nucl. Fusion **60** 126040 (2020).

[36] R.G. Kleva, J.F. Drake, "Density limit disruptions in tokamaks", Phys. Fluids B **3** 372 (1991).

[37] A. N. Kolmogorov, "On the Conservation of Conditionally Periodic Motions under Small Perturbation of the Hamiltonian" Dokl. Akad. Nauk SSR **98** (1954).

[38] J. Moser, "On invariant curves of area-preserving mappings of an annulus," Nachr. Akad. Wiss. Göttingen Math.-Phys. Kl. II 1962 (1962), 1-20.

[39] V. I. Arnold, "Proof of a theorem of A. N. Kolmogorov on the preservation of conditionally periodic motions under a small perturbation of the Hamiltonian" Russ. Math. Surv. **18**, 9-36 (1963)

[40] R.J. Hastie and T.C. Hender, "Toroidal internal kink stability in tokamaks with ultra flat q profiles", Nucl. Fusion **28** 585 (1988);

[41] F. J. Artola, K. Lackner, G. T. A. Huijsmans et al., "Understanding the reduction of the edge safety factor during hot VDEs and fast edge cooling events", Phys. Plasmas **27**, 032501 (2020);

[42] H. M. Smith and E. Verwichte, "Hot tail runaway electron generation in tokamak disruptions", Phys. Plasmas **15** 072502 (2008).

[43] P. Aleynikov and B. N. Breizman, "Generation of runaway electrons during the thermal quench in tokamaks", Nucl. Fusion **57** 046009 (2017).

[44] D. Hu and C. Liu, "Hot-tail electrons' impact on assimilation and injection penetration of D2 shattered pellet injections" Nucl. Fusion **62** 096002 (2022).

[45] H. W. Müller, K. Büchl, M. Kaufmann et al., "High-β Plasmoid Drift during Pellet Injection into Tokamaks" Phys. Rev. Lett. **83**, 2199 (1999).

[46] A. Matsuyama, R. Sakamoto, R. Yasuhara et al., "Enhanced Material Assimilation in a Toroidal Plasma Using Mixed $H_2 + Ne$ Pellet Injection and Implications to ITER" Phys. Rev. Lett. **129**, 255001 (2022).

[47] S. Futatani, G. Huijsmans, A. Loarte, et al., "Non-linear MHD modelling of ELM triggering by pellet injection in DIII-D and implications for ITER". Nucl. Fusion **54** 073008 (2014);

[48] T.C. Hender, J.C Wesley, J. Bialek et al., "ITER Progress Chapter 3: MHD stability, operational limits and disruptions" Nucl. Fusion **47** S128-S202 (2007).