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

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 magneto-hydrodynamic (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.

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 magneto-hydrodynamic (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. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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][2][3][4] and simulations [5][6][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][11][12][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 500 m s −1 injection velocity, 5% neon mixture ratio and reasonable fragmentation size, it takes about 10 ms 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 (100 eV), 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 ms 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.

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.

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 (CE) impurity model detailed in [6,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, ϕ) for the completeness of the paper, without going into great details: In the above equations, the Poisson bracket {f, g} ≡ R (∇f × ∇g) · ∇ϕ and ∇ pol represents the gradient within the poloidal plane. Equations (1) and (2) give the magnetic and velocity fields. Equation (5) is the induction equation, and η is the local Spitzer resistivity [21]. Note that the effective charge contribution [22] is included Here the threshold temperature T thres = 1 eV, below which the resistivity no longer varies as a function of the temperature. Further, equations (6) and (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 , ϕ f ): In this paper we choose ∆r NG = 0.08 m. Although the individual density source width from this choice may be underresolved 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.5 m, thus the ablation density source is generally well resolved. This also makes the result insensitive to the exact choice of ∆r NG so long as it is not larger than the plume width. Our toroidal resolution limit prevents us from using realistic ∆ϕ 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 ∆ϕ NG , we test two values: ∆ϕ 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 [23]. 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 [24]. 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: in which the parameter λ is a function of the mixture ratio X ≡ ND2 NNe+ND2 : r p is the fragment radius and ζ B ≃ ) 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 [23]. 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 [24,25]. 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 equations (8) and (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, equations (10) and (11) are the ion and electron pressure equations. The Braginskii parallel heat conduction [26] is used for each species respectively. With the initial density, the perpendicular particle diffusion is set to be constant at 6 m 2 s −1 and the perpendicular heat conduction at 10 m 2 s −1 . The perpendicular viscosity is set to 30 m 2 s −1 and the parallel one to 150 m 2 s −1 . The collisional thermalization between the ions and the electrons is represented by (∂ t T i ) c,e and (∂ t T e ) c,i [6]: where ν e/imp and ν 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 [27]. The charge state evolution is then determined by [19]: Here f i is the fractional abundance, S i and α 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 [28] P rad = n e n imp L rad (n e , T e ) ) .
Here n imp is the impurity density including all charge states while n i+ imp is that of each individual state, L i+ L , L i+ R and L i+ B 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 [6].

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.3 T and the total plasma current is I p = 15 MA, the central pressure is P 0 = 8.654 × 10 4 Pa, the total initial thermal energy is E th (0) = 36 MJ. The initial equilibrium profiles, including the electron temperature, density, the toroidal current density and the safety factor q profile, are shown in figure 1. This equilibrium exhibits a slightly reversed shear near the magnetic axis, and a significant portion of the core plasma (up to normalized fluxψ ≃ 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 figure 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 ∆ϕ NG = 0.3 cases and up to the third order in strength for the ∆ϕ NG = 0.5 cases. Ideal boundary condition is assumed in this study. We call this equilibrium the 'standard equilibrium'.

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 • , 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 = 200 m s −1 and the 'fast case' with v p,0 = 500 m s −1 . The fragment size distribution follows that of the Statistical Fragmentation model [33,34]: 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, κ 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 equation (18) over r p , we get: 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 × 10 22 neon atoms and 1.8 × 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 κ −1 p ≃ 1.43 mm. 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.

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 fragments 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 ∝ 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: κ −1 p ∝ L. Combined with equation (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 4/3 p , so that the total ablation rate from all fragments should have the following relationship with the injection amount: 0 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 ∝ R 3 is the plasma volume), then we have τ abl ∝ 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 [14], such effect could result in electron temperature collapse time on the order of 10 ms 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 2 η dependence in the resistive diffusion time where L η is the length scale of current redistribution. Meanwhile, the fragment flying time is only linearly depending on L η . 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.

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 sections 2.1 and 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.

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 500 m s −1 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 ∆ϕ NG = 0.5. We call this the 'reference case'. The SPI simulations with the same equilibrium but different injection velocities, ablation models and ∆ϕ 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 figure 3, covering the pre-TQ and early TQ phase. In agreement with previous CE impurity results [6], the temperature profile shows a 'scrapingoff' 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 figure 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 figure 4(a). For comparison, the MHD spectrum of a 200 m s −1 case with the 'non-local' ablation scheme as introduced in section 2.1 is also shown in figure 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 figure 4(b) simply because we run that case longer, 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 500 m s −1 case shows weaker MHD amplitude compared with the 200 m s −1 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 10 eV 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 sizeeffect we discussed, but takes several milliseconds instead, as is shown in figure 5 for the 200 m s −1 case with ∆ϕ 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 figure 4(b). We hereby define the beginning of the core relaxation as the beginning of the TQ, but we remind the readers that it does not corresponds 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 6 MJ of difference between the thermal energy loss and the radiated energy at t ≃ 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 1 MJ 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, 25 MJ out of the initial 36 MJ thermal energy is lost. The strong n = 1 MHD response is reflected in the density transport. The convective density transport at the time of the TQ onset is shown in figure 6. At the time of the TQ onset, significant density deposition exists on the q = 1 surface as marked by the density peak in figure 6(a). As the 1/1 kink develops, density cells can be seen 'dragging' into the q = 1 surface as shown in figures 6(b) and (c). The density cells were transported over a distance of O (1 m) in 0.2 ms, during which the fragments only travel 0.1 m. 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 500 m s −1 non-local ablation 3D case electron density and that of an identical simulation but in 2D (assuming axisymmetry) is presented in figure 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 in seen in the 3D case as is shown by the green curve comparison between figures 7(a) and (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 Kolmogorov-Arnold-Moser (KAM) surfaces [37][38][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 thus stronger drive for the instabilities.

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 [23].
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 equation (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 equation (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 200 m s −1 non-local ablation case are shown in figure 8. The different colours represent profiles at different times, while the solid lines and dashed lines represent the electron and the ion temperature, respectively. 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 eV) range as is shown in figure 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 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 m and R = 7.5 m can be seen at t = 2.12 ms and t = 2.45 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 figure 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 figure 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 is shown in figure 9(b), the 200 m s −1 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 figure 10, where the toroidal current density J ϕ profile multiplied by the major radius R is shown for both the 500 m s −1 (reference case) and the 200 m s −1 cases. Comparing figures 8(a) and 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      figure 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 500 m s −1 case, although the MHD amplitude is stronger for the 200 m s −1 case. As the TQ proceeds, the radiative collapse and the current redistribution further develop, the edge current becomes more perturbed.
Here one begin 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 figure 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 local 'electron temperature holes'. An example of such structure during the pre-TQ phase at t = 1.26 ms for the 500 m s −1 , non-local ablation case is demonstrated in figures 11 and 12.
In figure 11, a comparison of the log 10 T e between our Braginskii parallel thermal conduction case ( figure 11(a)) and a constant thermal conduction case ( figure 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 500 eV Braginskii value. The most notable feature in the Braginskii case is a dark blue spot around R ≃ 7.75 m and Z ≃ 0.65 m in figure 11(a), which is absent in figure 11(b). Such a 'electron temperature hole' corresponds to the 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 figure 11 is shown in figure 12, where the electron temperature ( figure 12(a)) and the electron density ( figure 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 figure 12(a), one can see the temperature hole develop 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 sometime 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 figure 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 /ν ee (n e , T e ) = ⟨σ⟩ −1 ≃ 2 × 10 21 m −2 .
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 figure 12(b) from δϕ = 0 to δϕ = 1, approximating the major radius as R ≃ 6 m. The resulting line integrated density is´n e dl ≃ 2.4 × 10 21 m −2 , comparable with the ⟨σ⟩ −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][43][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 [24,25]. 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 figures 11 and 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.

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 figures 3 and 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 500 m s −1 'reference' case studied in this section, and also for the  figure 13(a) is the 500 m s −1 , local ablation scheme case with ∆ϕ 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 figure 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 figures 3 and 8. This 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][38][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 figure 13(b), the same plot is shown for the 200 m s −1 , non-local ablation scheme case with ∆ϕ NG = 0.3 at t = 4.85 ms, shortly after the core relaxation. The edge field lines shown in figures 13(b) and (e) are more stochastic than those of figures 13(a) and (d), as is suggested by their better mixing and shorter connection length. This is due to the stronger MHD amplitude as is shown in figure 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 figure 13(b) and the persistent, although reduced, well-confined region in figure 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 figures 13(c) and (f ), representing the 200 m s −1 , non-local ablation scheme case at t = 6.36 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 m. Despite not running long enough in the 500 m s −1 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 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 m or R init = 7.5 m and Z init = 0.6 m for the 200 m s −1 case discussed above. The resultant Poincaré plot is shown in (θ, ψ  The core stochasticity is stronger since the blue curve reaches its steady result within 100 toroidal turns in both figures 15(a) and (b), while it takes hundreds of toroidal turns for the red curve to spread out as indicated by the slowly changing mean position and its variance towards their steady value. On the other hand, at t = 6.05 ms, the level of stochasticity is strong both in the core and the edge. The red and the blue curve quickly converge to the same mean position within several tens of toroidal turns as is shown in figure 15(c), and similarly for the variance in figure 15(d). Both steady values correspond to a spread over the whole volume, indicating good stochastic mixing. Such stochasticity evolution has important impact on the early-TQ heat flux dynamics as well as later seed runaway electron loss dynamics, which we will discuss in section 5.

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.

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 figure 16. All cases share the same velocity spread and the same 500 m s −1 reference injection velocity, and the two 3D cases share the same density toroidal deposition width ∆ϕ NG = 0.5. It can be seen from figure 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 figure 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.5 ms, 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 figure 16(b), so that only the density localization effect remains in the total radiation power, resulting in the stronger 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). radiation for the red curve shown in figure 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.75 ms and t = 2.1 ms in figure 16 roughly corresponds to the MHD activities at the same time in figure 4. We do not see significant difference in the ionization energy in general due to the robustness in the neon assimilation shown in figure 16(b).
Comparing figures 16(a) and (b), one may find it counterintuitive 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 energy and the total thermal energy loss, as is shown in figure 17. There, different colours represent the same cases as is shown in figure 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 nonlocal 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 figure 5, as time goes on this radiated power fraction increases.
We now move on to investigate the effect of the magnetic constraining term ζ B in equation (14). Two cases with ) 0.843 dependence. Such behaviour is consistent with our previous discussion concerning figure 16, that enhancement in the ablation rate would lead to only limited improvement in the assimilation due to its selfregulated nature, while the impact on the radiation could be more pronounced.

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 500 m s −1 , 3D non-local ablation case with ∆ϕ NG = 0.3 against the three cases discussed in figure 16 where ∆ϕ NG = 0.5. The result is shown in figure 19.
As expected, the enhanced toroidal localization indeed causes stronger radiation power as can be seen by comparing the magenta and red curves in figure 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 figure 19(b). Meanwhile, the radiation power increase is less than proportional to the change in the ∆ϕ 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.
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.

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 × 10 17 m −3 and tungsten density n W = 3.5 × 10 11 m −3 , according to the expectation of such L-mode discharge [29]. The background impurities are modelled using the CE assumption. Comparisons between the total impurity radiation power and that from the background species alone are shown in figure 20 for the 500 m s −1 case and the 200 m s −1 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% of the background impurity radiation. Although we consider a uniform background impurity distribution here, a nonuniform impurity distribution keeping 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.

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 figure 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 500 m s −1 , non-local ablation, ∆ϕ NG = 0.3 case with constant thermal conduction corresponding to that of the Braginskii value at T e = 500 eV in figure 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 figure 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 figure 22, which is the same as figure 19 but with an additional cyan curve that represents our constant thermal conduction case. Apart from the similar assimilation shown in figure 22(b), the constant thermal conduction case shows a delayed rise of the radiation power as can be seen in figure 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 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.

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 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 500 m s −1 fast injection case. The characteristic postinjection 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 CE 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 figures 4(b), 5 and 13, 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-tovolume 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 1000 m s −1 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 mixtureratio 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.