On the role of preexisting MHD activity for the plasma response to massive deuterium injection

Major disruptions are a serious threat to the successful operation of ITER because they lead to mechanical and thermal heat loads, the generation of runaway electron (RE) beams, and vertical displacement events (VDEs). A robust and thoroughly validated disruption mitigation scheme (DMS) in the case of failure is unavoidable and remains one of the outstanding challenges for ITER.¹ 

Massive material injection (MMI) is seen as the main strategy for disruption mitigation. In the past, the focus was mainly on massive gas injection (MGI). Due to its deeper plasma penetration and, hence, more uniform cooling, mitigation by shattered pellet injection (SPI) is now planned for ITER. In the case of deuterium injection, unshattered or weakly shattered pellets might also be a promising option for achieving material assimilation deeper inside the plasma.

An SPI system has a large set of degrees of freedom, such as the pellet material and size, injection schemes or position and numbers of injectors. Optimal values for some of these parameters will depend on plasma parameters as well as the focus of mitigation. Finding the optimal parameters for effective and, at the same time, highly reliable mitigation schemes is one of the subjects of ongoing work: one scheme, which targets the avoidance of runaway seeding, proposes diluting the plasma with deuterium SPI² in the early phase of the disruption as a first stage of a two-part mitigation. In this scheme, impurity SPI is applied only in a second step to mitigate thermal heat loads during the thermal quench (TQ). Numerical studies have shown that this fast dilution is possible without immediately triggering strong MHD activity, which would lead to a TQ. Here, a late TQ or even no TQ (or, in other words, a long cooling time) is favorable for a number of reasons: first, the cold post-TQ temperatures would massively inhibit the shard ablation making the dilution process less effective. Second, time is needed to trigger the impurity SPI in the second step. The third, and possibly the most important reason, for achieving a slow and gentle cooling of the plasma is that this is likely to and up with a cold Maxwellian population with no hot tail component, hence avoiding hot tail RE generation. A first set of predictive deuterium SPI simulations for an ASDEX Upgrade H-mode case with the variation of injection and plasma parameters has recently been conducted,³ demonstrating dominating modes with toroidal mode number n = 1 but also demonstrating the excitation of higher order modes and discussing cases in which a TQ was triggered by the injection of at least 10²¹ atoms. In those simulations, SPI was applied to a “healthy plasma” without MHD activity before the injection.

However, in practice, the mitigation system might be triggered when the plasma is already MHD active, e.g., an $m/n=2/1$ neoclassical tearing mode⁴ (with poloidal mode number m) may form a preexisting island structure. Simulations of impurity MGI with preexisting islands were performed previously for J-TEXT^5,6 and DIII-D.⁷ 

This paper presents the first global simulations for deuterium (shattered) pellet injections into an MHD active plasma carried out with the non-linear extended MHD code JOREK.^8,9 The focus of this work is on the aforementioned strategy of plasma dilution for RE suppression, where our aim is to understand whether it is still applicable in a pre-TQ phase. Injection parameters that are marginal for triggering a TQ shortly (within $∼1 ms$⁠) after plasma dilution has happened are of great interest. These cases can be seen as a worst-case, where the TQ occurs just after the plasma density has been increased. A possible worsening effect of preexisting islands is most critical in this scenario. Due to the impact of the islands on the local temperatures around the shard cloud, they may have a direct impact on the ablation. If it leads to too low temperatures in inner regions, full ablation not be ensured any more, which would lower the effectiveness of plasma dilution. In turn, higher initial perturbation amplitudes may increase the risk of an early TQ, when they superpose with the perturbations triggered by SPI. This could shorten the cooling time to values not sufficient for plasma dilution. Hence, the following central questions are addressed in this paper:

• What are the effects of preexisting islands on TQ formation?

• Is plasma dilution prior to the TQ also viable using deuterium injection in MHD active plasmas?

In this study, an ASDEX-Upgrade L-mode-like plasma is considered. As basic physical mechanisms of mode growth and the TQ triggering induced by massive material injection are investigated here and detailed experimental comparisons are left for further work, the more challenging dynamics with a steep H-mode pedestal are avoided. These mechanisms are studied by varying the initial island size w_i, number of injected atoms $NSPI$⁠, and the relative phase between island O-point and injection location $ϕO$⁠.

This article is structured as follows: Sec. II describes the simulation code, the simulation setup, and the simulation parameters used for this study. Section III presents the simulation results of various injection scenarios considered along with a detailed analysis. Finally, Sec. IV provides a concise summary and interpretation of the work, draws conclusions regarding implications for ITER, and provides an indication of future work needed to fully answer key questions for ITER.

The simulations performed in this study are carried out with the non-linear extended MHD code JOREK.^8,9 A reduced MHD model with fixed boundary conditions is applied, that includes extensions for a neutral deuterium fluid and (shattered) pellet ablation. A full description of the model is given in Refs. 9 and 10. For convenience, we repeat the equations here,

The eight simulation variables are the poloidal flux ψ, the poloidal flow potential u, the toroidal vorticity ω, the toroidal current density times the major radius j, the electron density ρ, the total temperature $T=Ti+Te$⁠, the parallel velocity $v∥$⁠, and the neutral mass density ρ_n. The resistivity is modeled according to the Spitzer¹¹ temperature dependency $η=η0(Te/Te,0)−3/2$ with the value on axis of $η0=1× 10−7 Ω m$ (a realistic value according to Spitzer, which also considers neoclassical corrections for electron trapping around the q = 2 surface, would be smaller by a factor of 2.5 according to Ref. 12, Chap. 14.¹⁰ Also, the parallel heat conductivity is temperature-dependent based on the Spitzer–Härm formula,¹¹ $χ∥=χ∥,0(Te/Te,0)5/2$ with a value on the axis of $χ∥,0=1.6×1011 m2/s$⁠,which is a realistic value within the experimental error of the electron temperature profile. The perpendicular heat conductivity is given as an numerical $ΨN$-dependent profile with values around $χ⊥=1.6 ± 0.3 m2/s$⁠. The dynamic viscosity is given by $μ=μ0(T/T0)−3/2$⁠, which is equal to $2.3×10−7 kg m−1s−1$⁠. The perpendicular particle diffusity is also given as an numerical profile with values around $D∥=2 ± 1 m2/s$⁠. The diffusion of neutrals is isotropic and has a higher spatially constant value $Dn=34 m2/s$ to mimic their ability to freely move independent of the magnetic field lines to some extent.

The poloidal plane is discretized by a flux-aligned 2D Bézier finite element grid. It comprises about 8200 elements. The toroidal direction is Fourier decomposed, and all modes from the axisymmetric part n = 0 up to $nmax=6$ are included. To investigate toroidal convergence, a case has been repeated with $nmax=10$ and is documented later in Sec. III A 1. Since the present article addresses the dynamics up to the TQ only, RE and free boundary extensions are not considered in this paper. For simplicity, background impurities are neglected, but their role could be important and this will be the subject of further studies. We also neglect Ohmic heating, which we consider because as a result of the absence of impurities, the plasma never becomes cold enough for Ohmic heating to make an important difference. Ohmic heating would be an important aspect if our aim was realistic modeling of flux surface healing. To reduce complexity, stabilizing background plasma flows are not considered—this corresponds to the assumption of an MHD mode being already locked to the vessel, this is the worst case because it is one of the most important causes of disruptions.^13,14 The bootstrap current source j_bs evolves for every time step based on the density and temperature evolution using the Sauter formula.^9,15 Thus, in a steady state limit, the Sauter bootstrap current density would be reached. This allows simulation of the neoclassical drive inside an island due to pressure flattening, although the Sauter formula is being used outside its original validity limits (axisymmetry, steady state). For simplicity, electron and ion temperatures are assumed to be equal, an assumption which becomes harder to justify during the TQ evolution that is not considered here.

Simulations are based on an ASDEX Upgrade L-mode-like equilibrium, as a pre-TQ MHD-active plasma has often already undergone H–L transition. The core electron temperature and density are at $Te,0=1.7 keV$ and $ne,0=2.1×1019 m−3$⁠, respectively. Temperature, current density, and q profiles are presented in Fig. 1, where the radial coordinate is the normalized poloidal flux defined by $ΨN=(Ψ−Ψaxis)/(Ψseparatrix−Ψaxis)$⁠. Below $ΨN$ always refers to the time varying flux surfaces.

The initial safety factor in the core, $q0=1.4$⁠, is taken to be well above unity, guaranteeing no intrinsic core instabilities. Some simulations with q₀ closer to 1 are performed for comparison and are summarized later in this article. Equilibrium is classically stable against tearing modes. The rational surface q = 2 is initially at $ΨN=0.59$⁠. A similar configuration is investigated in Ref. 16.

While most simulations could be continued through the TQ from a numerical point of view, we do not simulate the actual TQ in this study to save on computational costs, since the questions of interest for this study do not require the details of the TQ to be analyzed. We define the TQ here as the point at which the plasma becomes fully stochastic and the temperature throughout the plasma flattens out at a low value of $<100 eV$⁠. Note that for a fully realistic TQ modeling, Ohmic heating, background impurity radiation, and separate electron and ion temperatures would have to be taken into account.

To initialize preexisting magnetic islands of different sizes without having to reproduce the true experimental triggering mechanisms for seed islands, a separate simulation is performed. For this, the current gradient has been increased around the q = 2 rational surface to render the 2/1 tearing mode classically unstable. Running this setup in a non-linear simulation separately, leads, after a linear growth-phase, to a saturation of the 2/1 mode at a width of around $7 cm$⁠. The non-axisymmetric components of all eight variables were extracted at five different points in time prior to the saturation phase, each representing perturbations of a 2/1 island (and its sideband harmonics) at a different amplitude. These perturbations are imported into the original, classically stable equilibrium. Restarting these five configurations leads to some rearrangements, where the amplitudes of each island reduce about 10% within less than 1 ms and become virtually constant after that. By doing so, we eventually obtain different initial island sizes of $0.8, 1.5, 2.5, 4.5,$ and $6 cm$⁠, respectively. After initializing this perturbation, the island is evolved according to the extended MHD equations as a neoclassical tearing mode and becomes persistent for at least the first 10 ms, which is longer than the timescale that we are interested in. Further neoclassical effects other than the bootstrap current are not included here.

For this kind of mitigation scheme, we consider a configuration of a “weakly shattered” pellet that is broken into just ten shards of equal size, allowing for deep penetration of the material. A strongly shattered deuterium pellet would instead produce a cloud of extremely small pieces and many gases, leading to much worst penetration properties and possibly no efficient core dilution. The pellet nozzle is positioned on the outer midplane and activated at $t=0.1 ms$—shortly after importing the magnetic island perturbation into the simulation.

The ablation is modeled by neutral gas shielding,¹⁷ and the neutral gas source term of every single shard is assumed to be of Gaussian shape both in poloidal and toroidal directions. The implementation is the same as in Ref. 10. In the present study, the poloidal radius of each neutral cloud is set to $8 cm$ for numerical reasons; however, this is unrealistically large. Generally speaking, a cloud radius of a few centimeters is reasonable: in the parallel direction, this is caused by the fast so-called ambipolar expansion direction, driven by the strong local electron pressure gradient and the fact that the plasmoid is transparent to the ambient field lines.¹⁸ In the perpendicular direction, this expansion already stops as soon as a radius of about $1 cm$⁠. However, a strong drift of the cloud toward the low field side can again produce stretching of a few centimeters in the perpendicular direction.¹⁹ The initial velocity vector of each shard is varied randomly both in magnitude and the direction around a reference velocity of $250 m/s$ pointing in the major-radial direction with a maximum variance from the reference speed of $40 m/s$⁠. The directions are chosen randomly such that all shards travel in a cone with aperture $α=0.349 radians$ around the reference velocity vector. These shard velocity vectors are identical in all simulations to allow for direct comparison. The averaged poloidal position of the shards at some time points are given in Fig. 1. The whole neutral source generated by all ten shards altogether is one large single cloud with a poloidal expansion of around $ΔR×ΔZ=20×16 cm2$ and toroidal expansion of $Δϕ=0.7$ [at $t=1 ms$⁠, measured at the contour line, where $ρn=max(ρn)/e$ with e being Euler's number]. Therefore, the shard cloud is always larger than the imposed island structures.

The relative phase of the injection is chosen in most simulations such that the shards directly hit the island O-point, which we label as relative toroidal phase $ϕO=0°$⁠. Variations of the phase are presented later on in this paper to assess the effect of different phases, where $ϕO=180°$ corresponds to X-point injection.

For the amount of injected atoms $NSPI$⁠, four different values are considered: $3×1020,6×1020,1×1021$ and $6×1021$⁠. For comparison, the initial total plasma particle content is $1.9×1020$⁠. Altogether, including simulations without a preexisting island and a modified q-profile, more than 70 non-linear simulations are performed and analyzed. One simulation is usually executed on more than 300 cores and takes up to 240 wall-clock hours to complete. This leads to a consumption of more than 2 Mcore-hours in total.

During the following analysis, we need to distinguish between shards only by diluting and cooling the plasma adiabatically and heat losses by degraded confinement. This is done by comparing the evolution of the total thermal energy $Etotal$ content as well as the total particle content $Ntotal$⁠. In a pure-dilution case, $Ntotal$ increases, while $Etotal$ remains constant. In turn, a degradation of the energy confinement results in a decrease in $Etotal$⁠. To improve the understanding of the local dynamics in the vicinity of the magnetic axis, we analyze the quantities “core particle content” (⁠$Ncore$⁠) and “core thermal energy” (⁠$Ecore$⁠). The core region is defined as the volume inside the original $ΨN=0.1$ surface and remains constant over time.

The evolution of island widths is estimated based on the formula

where c = 0.7 is a prefactor determined empirically by comparing the analytical expression with Poincaré plots at several time points. $Ψm,n$ denotes the m/n component of the poloidal magnetic flux in straight field line coordinates on the rational surface. A similar approach has been applied, e.g., in Ref. 20. This makes it possible to “measure” the island width continuously throughout a simulation and can continue to be applied in stochastic scenarios in which the Poincaré plots no longer allow measurements of the island sizes.

In the following, we describe in sequence the simulations performed to investigate the interplay of preexisting MHD modes with massive deuterium injection. First, in Sec. III A, we will focus on injections with $NSPI=6×1020$⁠. We find that this amount of material is needed to trigger a TQ shortly after core dilution, once all the material has been ablated. As indicated in Sec. I, this scenario is one of the primary motivators behind this study. It can be considered a limiting case—a threshold—between the two regimes in which there is an immediate TQ before full ablation or in which no TQ occurs after injection. We expect the influence of preexisting MHD activity to be most critical in a scenario close to this threshold. The case with $NSPI=6×1020$ injected into an unperturbed plasma is studied is detail in Sec. III A 1 to establish a baseline for subsequent comparisons. Injections into the island O-point of preexisting magnetic islands of different initial island sizes are studied in Sec. III A 2. We see here that for the case of a large preexisting island, the plasma response is significantly affected. The dependency on the injection phase with respect to the preexisting island is studied in detail in Sec. III A 3, where the injection location relative to the island is systematically varied. Here, we highlight, in particular, that an injection into the immediate vicinity of the island X-point leads to an earlier TQ than the injection into other phases. A preliminary not yet fully developed explanation for different TQ timings is presented in Sec. III A 4. Section III B deals with variations of the amount of injected material. We first discuss the effect of that variation for the injection into an unperturbed plasma in Sec. III B 1. We analyze the effect of different amounts of material injected into the O-point or X-point of a large preexisting island in Sec. III B 2. These results suggest that the injection close to the X-point seems to lower the threshold of the amount of material injected to trigger a TQ. To study how sensitive our results are to the q-profile, injections into a plasma with a modified q-profile are discussed very briefly in Sec. III C.

The plasma equilibrium used for our studies initially contains a thermal energy of $Etotal,0=71 kJ$ and $Ntotal,0=1.9×1020$ deuterium ions. The core region (defined by $ΨN≤0.1$⁠) initially contains a thermal energy of $Ecore,0=14.5 kJ$ and $Ncore,0=0.2×1020$ deuterium ions. In the absence of an injection by diffusive cross field transport, 9.5% of the thermal energy and 0.5% of the total particle content are lost within 4 ms (see Fig. 2). For simplicity, no other sources are involved other than the neutral source from the pellet injection, which could compensate for the above losses.

An overview of the dynamics with $NSPI=6×1020$ is given in Fig. 2 where the time evolution of the particle content and thermal energy content in the complete plasma and in the core region are shown. The shards are injected at $t=0.1 ms$ and begin to affect the plasma dynamics at a simulation time of $t=0.3 ms$ by increasing the total particle content as well as by decreasing the total thermal energy content. Full ablation is reached only at $t=1.9 ms$⁠, and the particle content increases by nearly 300% with a high fraction of material assimilated around 93%. In the first millisecond, $Etotal$ decreases by around 10% due to losses induced by edge MHD activity and then decays only slightly, indicating a restored edge confinement. In comparison to previous deuterium-SPI studies for AUG,³ the material assimilation is strongly enhanced because our study concerns an L-mode plasma that is less prone to violent ELM-like edge instabilities. A rapid drop of $Ecore$ by more than 50% takes place within about $100 μs$ at around $t=1.2 ms$⁠, while $Ncore$ increases less abruptly due to different time scales for the radial transport of heat and particles. After this core thermal energy drop, the total and core thermal energies remain almost constant for $1.4 ms$⁠, while $Ncore$ increases monotonically, reaching a maximum value of about $0.8×1020$⁠, which is an increase in 400%. After $t=2.5 ms$⁠, total and core thermal energy, as well as core particle content drop, which marks the beginning of the TQ as indicated by full stochastization of the plasma. We remind that we do not continue the simulations beyond this point. As this case, involving a delayed TQ shortly after plasma dilution, is of particular interest in this study, it will be analyzed in detail in the rest of this subsection.

The time evolution of the magnetic perturbation spectrum excited by the material injection is shown in Fig. 3. When the shards penetrate the plasma, all toroidal $n≠0$ modes included in the simulation (⁠$n=1,…,6$⁠) are excited, as expected for MMI-simulations. The n = 1 mode is dominant throughout most of the simulation because an injection at a single toroidal location is under investigation. At around $t=0.8 ms$ when the shards arrive at the q = 2 surface, the n = 1 magnetic perturbation is transiently decreasing and the n = 2 mode becomes dominant for a few hundred microseconds. The reason for this monotonic evolution of the harmonics is associated with a phase transition of the 2/1 perturbation and will be discussed later in this section. When $Ncore$ begins to increase after $t=1 ms$⁠, we start to see a strong excitation of the fifth harmonic starting from low amplitude. Slightly later, the other n > 1 modes also begin to grow quickly due to strong non-linear mode coupling, which eventually leads to a burst of MHD activity at around $t=1.2 ms$ related to the first crash that decreases core confinement. At that time, the shards are located around the q = 1.5 surface. During this crash, there is a strong increase in $Ncore$ as well as a partial loss of $Ecore$⁠, as described above. However, the total particle and thermal energy content of the plasma as well as the ablation rate are only weakly affected, indicating that this crash corresponds largely to a re-distribution inside the plasma.

Subsequently, the core is successively further diluted while all harmonics except for n = 1 drop transiently in amplitude. The shards are fully ablated at $t=1.9 ms$⁠, which implies that some shards reach the $ΨN=0.1$ surface. Later, at around $t=2.3 ms$⁠, an increase in modes of higher order becomes visible again and a distinct mode coupling in the next several hundred microseconds marks the second crash at around $t=2.5 ms$⁠, which is more violent than the first. During this crash, $Ncore$ and $Ecore$ decrease simultaneously and successively; from $t=2.6 ms$ on, the total energy starts to drop, indicating the TQ. The simulation is not continued further because the TQ dynamics are not part of the scope of this article.

The large n = 5 mode, found particularly during the first core crash, raises the question of whether convergence occurs in the toroidal resolution for this set of simulations. To check this, the same simulation was repeated with $nmax=10$⁠. The timing and the quality of both crashes remain similar in terms of the thermal energy and particle content evolution. Also, the overall evolution of the n = 1 magnetic energy does not differ significantly, indicating a similar behavior of the prominent 2/1 mode. During the first crash, however, n > 4 modes are still excited but do not attain such large amplitudes; thus, the first and second mode continue to dominate. As a result, perpendicular heat transport in the core is not as efficient and $Ecore$ stabilizes at a slightly higher level after the first crash. The behavior of all modes after the first crash is ultimately fairly similar. As higher modes only induce small changes in the overall dynamics, we conclude that the toroidal resolution of $nmax=6$ is sufficient for this study in which the general dynamics are primarily determined by low-order modes. With the help of this run, it is also possible to check the convergence of the modeling of the gas cloud. The size and shape of the cloud is virtually the same for both toroidal resolutions providing confidence that $nmax=6$ is sufficient for the parameter scans.

The dynamics described in the previous paragraphs are also reflected in the magnetic topology as shown in a series of Poincaré plots in Fig. 4. Just before the first crash (at $t=1.2 ms$⁠), the 2/1 island reaches a width of $w2/1=60 mm$ and, in particular, 3/2 and 4/3 modes are excited. During this first burst of MHD activity (⁠$t=1.2,…,1.55 ms$⁠), nearly the full region inside q = 2 becomes stochastic, while flux surfaces at the edge of the plasma remain intact. The 2/1 island has grown to $w2/1≈70 mm$ by the end of the first crash. This is followed by a reformation of flux surfaces in the center region for $ΨN<0.3$ (⁠$t=1.7,…,2.2 ms$⁠). Successively, also the outer regions become more stochastic, as the 2/1 and 3/1 modes continue to increase. Since the connection length of field lines from inside the plasma to the targets nevertheless remains long, which can be seen from the Poincaré plots (e.g., at $t=1.9 ms$⁠), $Etot$ continues to decrease slowly only. At the onset of the second crash, the width of the 2/1 island corresponds approximately to $w2/1≈170 mm$⁠. Finally, inner flux surfaces start to break up again and from $t>2.5 ms$⁠, the whole plasma becomes fully stochastic marking the TQ. Short connection lengths in this phase are reflected in a lower density of points in the Poincaré plot and lead to a faster decrease in $Etot$⁠.

The island phases with respect to the injection location are worth discussing. Before the shards pass the q = 2 surface, 2/1, 3/1, 4/1, and 5/1 islands arise with X-points in phase with the injection location. When the shards pass the q = 2 surface, the 2/1 island exhibits a phase shift of $Δθ=−90°$ in the poloidal direction and continues to grow, now with the O-point in phase with the injection location. The origin of this behavior is discussed below in this subsection.

In Fig. 5, the time evolution of the (toroidally averaged) density, temperature, and pressure profiles across the outer midplane is shown. Up to the first core crash, the effect of the shards on the temperature profile is primarily adiabatic cooling of the outer regions up to the q = 1.5 surface, while the pressure profile remains largely unchanged. Also, the central temperature and density remain constant. Consequently, a hollow density profile establishes up to this point in time. During the first crash, around $t=1.2 ms$⁠, the central temperature drops from $1400 eV$ to below $250 eV$ within $300 μs$ and becomes fairly uniform across the domain. This is mainly caused by enhanced radial heat transport along field lines in the stochastic core. As particle transport along field lines takes place over the far longer timescales of the ion sound speed, the density profile changes less quickly and remains hollow such that a hollow pressure profile also becomes established. As flux surfaces re-form after the first crash, radial transport reduces and the structure of a hollow pressure profile remains up to the second crash. The shards, which are now close to the magnetic axis, provide a continuous core fueling up to $t=2.1 ms$⁠, which leads to further adiabatic cooling. However, this effect is limited because of the small remaining shard sizes and reduced ablation due to low temperature. The outer region of the pressure profile relaxes (⁠$t=1.8,…,2.4 ms$⁠), as a result of increasing field line stochasticity. In the second crash, the entire pressure profile flattens (⁠$t=2.5,…,2.7 ms$⁠). The effect on the temperature profile is moderate because the temperature had already flattened.

At the beginning of the first crash, there is still only one identifiable single neutral cloud with a radial width of about $18 cm$⁠. As a result of the shard velocity distribution, the shards can be divided into an inner and outer subgroup, which are separated by a distance of $5 cm$⁠. As the inner shards face hotter electron temperatures, they ablate faster than the outer shards and the whole neutral cloud has a slightly larger particle density in the inner half. After $t=1.5 ms$⁠, only the “outer” shards remain and penetrate the cooled core much deeper, where they too fully ablate. A more uniform velocity distribution would probably lead to a more peaked neutral cloud and, possibly, a slightly faster core crash or less efficient core fueling after the crash. However, a change in the overall dynamics due to a different velocity distribution is not expected.

As discussed in previous publications,^10,21 there are two mechanisms of mode-destabilization by material injection, which are related to the temperature dependency of the Spitzer-like resistivity $η=η0(Te/Te,0)−3/2$⁠. The first mechanism is driven by local cooling on rational surfaces, which directly induces helical current perturbations via the helical perturbation of the resistivity. The second mechanisms is caused by an axisymmetric cooling, which increases η at the outer parts of the plasma such that the current profile contracts and possibly destabilizes modes via the classical tearing instability. To assess by which of the two mechanisms, the MHD modes are primarily driven, we performed a simulation in which η is only calculated from the axisymmetric n = 0 component of the temperature $Te,axi$⁠. This removes the resistivity perturbation from helical cooling and, hence, the related mechanism of mode-destabilization. Starting the simulations at $t=0.1 ms$ shows at first a similar behavior compared to the original case in terms of the evolution of the island widths (see Fig. 6). After $t=0.6 ms$⁠, the evolution bifurcates and the island growth stagnates at a width of $w2/1=18 mm$⁠, while no islands with higher n establish. This suggests that the helical cooling mechanism dominates the dynamics leading to the TQ in our simulations.

As pointed out above, in the early phase of injection, the shards trigger n = 1 islands of small width and X-points in phase with the injection location. We analyze this behavior based on the 2/1 mode below, where it is found to be the most pronounced. Similar observations have been made for islands of higher poloidal mode numbers. In this first phase, the $Ψ2/1$ component of the poloidal flux $Ψ$ exhibits a phase jump of π at the q = 2 surface (Fig. 7). The amplitude is close to zero on the rational surface (corresponding to a very small island size) and shows local maxima slightly inside or slightly outside the surface, which indicates a kink-parity dominated response. Only after $t=0.8 ms$⁠, when the shards have reached the q = 2 surface, the tearing structure begins to grow, interferes with the kink structure, and becomes dominant after around $t=1 ms$⁠. At this time, the phase jump disappears. Now, a 2/1 island with the O-point in phase with the shards becomes dominant (as expected for a magnetic island driven by helical cooling on the resonant surface). As expected, its width increases together with the amplitude of $Ψ2/1$ at the q = 2 surface.

The origin of the island phase prior to helical cooling is not fully clear but it is probably a response to a strong deformation of the equilibrium. $0.6 ms$ after injection, the ablation has led to a strong pressure perturbation, in particular, around the q = 2 surface. This results in a radial pressure gradient of $|∇pr|≈180 kN/m3$ on the low field side, which is a quintupling compared to the initial values (see Fig. 8). This produces a perturbed flux $Ψ2/1$ of kink parity as a response, leading to flux surface deformation and island formation via reconnection. Scanning through $NSPI$ shows further that $Ψ2/1$ scales proportionally to the amplitude of $∇pr$ around q = 2 at $t=0.5 ms$⁠. Up to this point, helical cooling is highly asymmetric and less pronounced on the high field side. As the process continues, the pressure gradient decreases as the material distributes and the shards move further inwards. After $t=0.9 ms$⁠, the temperature perturbation $Te∗$ then dominates and produces a 2/1 helical current perturbation, which drives the tearing mode.

To conclude this section, we propose that mode excitation is an interplay of an excited pressure perturbation due to a local increase in density and the resistivity perturbation due to helical cooling. In the first phase, where the pressure perturbation is very peaked, a kink mode is excited by strong pressure gradients, which leads to reconnection in the region of the shards. Only after some relaxation of the pressure profile and the onset of localized cooling at the q = 2 surface is a tearing mode is excited that produces a magnetic island with the O-point around the shards.

In the following, simulations will be repeated with the same setup and parameters but including preexisting 2/1 islands of different widths w_i. Section III A 2 investigates cases of injections into the O-point region of the island. Then simulations with different injection phases are analyzed in Sec. III A 3.

Having studied injections without preexisting islands in Sec. III A 1, we now turn to injections with $NSPI=6×1020$ into the O-point of a preexisting 2/1 NTM. A comparison in terms of particle and energy content of selected simulations has already been given in Fig. 2. The impact of the island on the ablation becomes visible for the first time in the evolution of $Ntot$⁠, when shards penetrate the O-point after $t=0.6 ms$⁠. Due to temperature flattening in the O-point, the ablation is enhanced outside the q = 2 surface and decreased inside. This faster or slower increase in $Ntot$ also has an impact on the excitation of edge instabilities, which is reflected by a different evolution of $Etot$⁠. Apart from this relatively weak effect, the evolution of $Ntot$ and $Ncore$ initially remains similar to the cases without a preexisting island, which we will now refer to as “no preexisting island cases” (abbreviated as “NPI-cases”).

The onset of the first core crash is only delayed by $0.1 ms$ for the largest island. Again, $Ecore$ becomes constant after $t=1.5 ms$⁠, but at a slightly larger value of $8 kJ$ independently of w_i. The preexisting island, however, has a strong influence onto the timing of full stochastization, i.e., the TQ. It occurs eventually with a delay of $1.1 ms$ at $t=3.6 ms$ for the $wi=6 cm$ case. As a consequence, both $Ecore$ and $Etotal$ remain constant for a significantly longer time compared to the NPI-cases. By contrast, the effect of the smaller island (⁠$wi=2.4 cm$⁠) on the TQ timing is very limited.

As the delay of the TQ for O-point injection with a large preexisting island is an interesting and potentially important observation, we now consider several different initial island sizes in more detail. The evolution of the island sizes in these cases is shown in Fig. 9. For all cases, the island width remains essentially constant up to $t=0.6 ms$⁠. After a short interplay between the preexisting island and the pressure-driven kink mode, which, in some cases, causes a slight shrinking of the island, all islands begin to grow distinctly after $t=0.8 ms$⁠. For larger w_i, the growth rates are smaller such that very similar island widths are observed at $t=1.2 ms$ for $wi≥4.5 cm$⁠. Shortly after this point in time, the first crash occurs in all cases, which begins with a strong growth of the n = 5 harmonic and then subsequently a non-linear coupling of higher harmonics around $t=1.2 ms$ (see Fig. 10). The evolution of these higher harmonics is fairly independent from the initial island width up to this point in time. The islands of smaller initial width $wi≤2.4 cm$ continue to increase monotonically with similar growth rates compared to the NPI-case. For all of these, the TQ begins at about $t=(2.6 ± 0.2) ms$ and the plasma again exhibits strong mode coupling [see Fig. 10(a)]. This is only $≈0.1 ms$ later compared to the NPI-case. In contrast, the behavior of the larger islands (⁠$wi=4.5 cm$ and $wi=6 cm$⁠) differs significantly from the cases with no or small initial island size. In these cases, the island sizes, as well as the n = 1 energy [see Fig. 10(b)], decrease transiently in the time window $1.2 ms<t<1.8 ms$⁠, which means that they even decrease during the first burst of MHD activity in the core. After that, they start to grow again. As a consequence of the decrease itself and the associated delay, a second crash sets in at just $t=3.2 ms$⁠. The mode coupling in the $wi=6 cm$ case is fairly weak compared to the NPI-case, and especially, modes of higher n remain at much smaller amplitudes such that this second crash does not lead to full stochastization. The TQ only occur at $t=3.6 ms$ as a third crash. The full stochastization corresponding to the TQ onset is reached at an island size of $≥17 cm$ in all cases. It should be noted that the island widths calculated at these time points can no longer be seen from Poincaré plots due to strong stochastization. As they remain a very intuitive and useful quantity for expressing the 2/1 perturbation amplitude, they are calculated instead from the 2/1 component of the poloidal magnetic flux at the q = 2 surface according to Eq. (9).

The decay of the first harmonic after the first crash, which we see in $Emag,1$⁠, as well as in $w2/1$⁠, is the main reason for a delayed full TQ. This will be investigated later in Sec. III A 3, where we compare this behavior to injections with other phases with respect to the island.

After a first look at the plasma dynamics with and without preexisting islands, we turn now to an analysis of the details. For this, we look, in particular, at the magnetic topology evolution by studying Poincaré plots for the case of the smaller preexisting island with $wi=2.4 cm$ (Fig. 11) and for the larger island with $wi=6 cm$ (Fig. 12). Although the smaller preexisting island does not influence the long-term dynamics after the first crash, dynamics in the early phase (⁠$t<1 ms$⁠) may be affected significantly: at around $t=0.5,…,0.8 ms$⁠, when the shards are close to the q = 2 surface, the initial island breaks up into a 4/2 structure implying that there is no transient 2/1 perturbation with tearing parity present. One X-point is then in phase with the shards, similarly to the X-point of the 2/1 mode in the early phase of NPI-cases. When the shards are inside the q = 2 surface (⁠$t=1.2 ms$⁠), the 2/1 structure has reformed with the same phase position as before, i.e., the island O-point is in phase with the shard location. A similar observation of such a transient 4/2 island was made in Ref. 7, where neon MGI was applied to a plasma with a preexisting 2/1 island. The further evolution is then similar to the NPI-case: after the first crash, flux surfaces in the center re-form and remain intact up to the second crash, while the outer regions continue to become more stochastic (⁠$t=1.7 ms$⁠). Full stochastization is reached at around $t=2.6 ms$⁠, hence not significantly later than in the NPI-case. We see, that this small island of $wi=2.4 cm$ does indeed have an impact on the magnetic topology up to the first crash, but the further evolution of the plasma is virtually the same as in the NPI-case. The transient disappearance of the 2/1 island that makes the underlying 4/2 structure visible is a result of the initial island with the O-point at the injection position canceling with the shard-induced island that has the X-point aligned with the shards. When the shards reach the q = 2 rational surface, helical cooling becomes the dominant drive of the island, which leads to regaining of the 2/1 structure with the O-point in phase with the shards (compare with Sec. III A 1).

In contrast, the overall plasma dynamics are strongly influenced in the case of a large preexisting island (⁠$wi=6 cm$⁠): when the shards are close to q = 2, a clear mode structure is not visible in the Poincaré plot as the surface becomes stochastic (⁠$t=0.6 ms$⁠). After that, the island re-forms and grows until $t=1.2 ms$⁠. The evolution inside the q = 2 surface is similar when compared to the other cases. During the first crash, the core becomes nearly fully stochastic (⁠$t=1.5 ms$⁠). This is consistent with the fact that the evolution of the magnetic energies of higher harmonics is similar in that time period for all w_i. After this first crash, the flux surfaces re-form in the inner region $ΨN<0.3$⁠, while the 2/1 island shrinks to about the initial width $w2/1≈6 cm$ at $t=1.8 ms$⁠. From this point on, the degree of stochastization increases and at $t=3.2 ms$⁠, just before the second crash, only remnants of the 2/1 and 3/1 islands are visible in the Poincaré plots. However, in the second crash, not all flux surfaces in the center become immediately stochastic: those for $ΨN<0.1$ still remain intact for a few hundred μs and vanish only at $t=3.6 ms$⁠, which corresponds to the full stochastization. The separation of the final phase before the TQ into two distinct crashes is only observed for the largest preexisting island. Together with the transient shrinking of the island, this delays the time of full stochastization (TQ) by more than $1 ms$ compared to the cases with no or only a small preexisting island.

After looking into the magnetic topology via Poincaré plots, we turn now to the evolution of the density, temperature, and pressure profiles for the large preexisting island of $wi=6 cm$ as shown in Fig. 13, which can be compared to the NPI-case shown in Fig. 5. Until the first crash, the evolution of temperature and pressure is very similar compared to the NPI-case with the exception of the behavior around the q = 2 surface: initially, a pronounced flattening of the temperature and pressure profiles caused by the 2/1 island located at $R=1,25 m$ and $R=2,01 m$ on the midplane can be seen. At $t=0.6 ms$⁠, this flattening degrades caused by the temporary stochastization of the q = 2 surface. Hence, the situation just before the first crash (⁠$t=1.2 ms$⁠) is nearly the same across all cases. The temperature in the center drops rapidly and a hollow pressure profile establishes during the crash. Until $t=1.8 ms$⁠, the electron temperature inside the q = 2 surface has flattens around $250 eV$⁠. This means that in the presence of the large preexisting island, the dynamics of the first crash take more time (about $300 μs$⁠) but are not drastically different from other cases. Subsequently, the evolution bifurcates: while for no or small islands, the hollow pressure profile starts to relax from $t=1.8 ms$ and relaxes completely within $0.7 ms$⁠, ending in the TQ; it remains hollow for a longer period for large w_i. It is only at around $t=2.3 ms$ that the start of pressure profile flattening can be seen. This delay of $0.5 ms$ is again related to the retarded evolution of the 2/1 and other n = 1 modes, which drive the stochastization of the outer regions. As a result of the second crash, the density becomes nearly completely flattened after $t=3.3 ms$⁠. However, a small maximum inside the core region (⁠$ΨN≤0.1$⁠) remains, and this only vanishes at full stochastization (⁠$t=3.6 ms$⁠).

We close this section with a discussion of the behavior of the 2/1 island just before the shards enter the q = 2 surface: at this early time point, we see distinctly different dynamics including a 2/1 mode with the X-point in phase with the shards (NPI-case), a 4/2 mode (case with small preexisting island), or a largely stochastic layer around q = 2 (case with large preexisting island). Which of these different phenomena are observed in a case, depends on the amplitudes of the initial perturbation $Ψ2/1,0$ and the perturbation $Ψ2/1,s$ that would be triggered by the shards alone and which is of opposite phase around q = 2. $Ψ2/1,s$ is identifiable from NPI-cases and is given for $NSPI=6×1020$ in Fig. 7. For cases with a preexisting island, we assume that the actual perturbation is given by a superposition in the early phase, $Ψ2/1(t)≈Ψ2/1,0+Ψ2/1,s(t)$⁠.

In this phase of the simulation, $Ψ2/1,s$ primarily represents a 2/1 kink with a small tearing component as discussed above. Its amplitude (at $ΨN=0.65$⁠) is $≈2.5×10−5 Wb$ at $t=0.6 ms$ with local phase $ϕ2/1≈0$⁠. The initial $Ψ2/1,0$ of $wi=2.4 cm$ has a similar amplitude at the same radial position (see Fig. 14, $t=0.2 ms$⁠) without a phase jump. As a consequence, both perturbations cancel at $t≈0.6 ms$ on the q = 2 surface and the otherwise subordinate $Ψ4/2$ component becomes visible as a 4/2 island structure (see Poincaré plots of Fig. 11). For the largest island, $Ψ2/1,0≈1.3×10−4 Wb$ at $ΨN=0.65$ is five times larger compared to $Ψ2/1,s(0.6 ms)$⁠. As a consequence, the initial tearing mode structure is only slightly impaired by the pressure perturbation and the small tearing component associated with the kink response. A cancelation is not possible in this case. During the later evolution in both cases, the tearing mode structures reestablish leading to the reformation of a 2/1 island with the O-point in phase with the shards. The dynamics in this phase are dominated by helical cooling.

Section III A 3 describes cases with injections at different phases with respect to the preexisting island. Then, the mechanisms by which large preexisting islands can affect the timing of the full TQ are discussed.

In agreement with the observations for injection into the O-point, which were discussed in Sec. III A 2 and are henceforth denoted “O-point cases,” the overall evolution of the total and core particle content prior to the TQ is not strongly affected by the preexisting island for other injection phases and is, therefore, not shown here. Nevertheless, an impact of preexisting islands on the energy quantities can be seen. An overview for the injection directly into the X-point (“X-point cases”) at different values of w_i is given in Fig. 15. The first (core) crash happens at exactly the same time $t=1.2 ms$ independent of w_i. However, the larger the initial island, the earlier the TQ occurs: for $wi=6 cm$⁠, this is the case at around $t=1.8 ms$⁠, which is reflected as a decrease in $Etot$⁠. This implies that the time between the first crash and the TQ has shortened to $0.6 ms$ and the TQ happens at a point in time, where shards are not fully ablated. $Ecore$ reduces to $6 kJ$ during the first crash and continues do drop only after $t=2.8 ms$⁠. This is a significantly lower value than in any other $NSPI=6×1020$ case in which $Ecore$ dropped to values between 7 and $8 kJ$ during the first crash. The following paragraphs describe these dynamics for the injection into the X-point of the large preexisting island in more detail.

The origin of the early TQ at $2 ms$ in the X-point case lies in the simultaneous strong destabilization of modes of higher order in the core as well the strong growth of n = 1 modes, as depicted in a series of Poincaré plots (Fig. 16). As soon as $t=1.2 ms$⁠, the outer region with q > 2 is much more stochastic compared to the situation in the O-point case (see Fig. 12). Together with the excitation of n = 2, 3 modes inside q = 2, the plasma is almost fully stochastic by $t=1.5 ms$⁠. Remnants of the 2/1 island disappear within the next $0.3 ms$⁠. That means, for the X-point cases, the first core crash turns into a TQ with full stochastization. In contrast, islands in the core have already reappeared at $t=1.8 ms$ in the O-point case.

The pronounced stochastization also affects the temperature and density profiles shown in Fig. 17:

• During the first crash (⁠$t=1.2,…,1.5 ms$⁠), the temperature becomes flattened over a wider radial range including the q = 2 surface and not just within the core. Consequently, the flattened electron temperature is $180 eV$ and significantly lower compared to the O-point or NPI-case.

• As the core remains fully stochastic afterwards (around $t=1.8 ms$⁠), stochastic transport is more efficient and a monotonic density profile is already established by then.

The result of both processes is that the hollow pressure profile decays already after the first crash and $Ecore$ remains at a smaller value.

Now that we have seen that an injection into the island X-point behaves very different from an injection into the O-point, we analyze the dynamics with different phases (between X- and O-point) below. The evolution of the island width for all phases under consideration is shown in Fig. 18. The dynamics in the X-point case (⁠$ϕO=180°$⁠) are enhanced; hence, a rapid island growth from $t=0.85 ms$ to $t=2 ms$ occurs. This is leading to a saturated island width of $>20 cm$ and eventually the earlier TQ identified previously. In the first phase, however, when the shards are on or close to the q = 2 surface (⁠$t=0.6,…,0.8 ms$⁠), the island width decreases transiently for a very short time. This is in contrast to the dynamics for the O-point case (⁠$ϕO=0°$⁠), where immediate island growth occurs when shards arrive at the q = 2 surface. From which two main phases are observed: in the first phase (⁠$t=0.6,…,1.0 ms$⁠), the island in the O-point case is larger than in the X-point case. Then, after a turning point at $t=1.0 ms$⁠, the situation reverses and the island in the O-point case is now smaller. Injecting between O- and X-points (⁠$0°<ϕO<180°$⁠) leads to a consistent behavior: the closer $ϕO$ gets to $180°$⁠, the greater the weakening in the first phase and the earlier the TQ. For $ϕO=90°$⁠, we also see a decrease in the size of the island during and shortly after the first crash (⁠$t≈1.1,…,1.7 ms$⁠), just as in the O-point case, which is, however, less pronounced. This decrease is not observed for values of $ϕO$ close to $180°$ (X-point phase).

In this subsection, we discuss the origin of the different behavior of the mode growth in the O-point, X-point, and NPI-cases, in particular, with respect to the different TQ times. The initial equilibrium is tearing mode stable mode and, as discussed on the basis of Fig. 6, the island growth prior to the first crash at $t=1.2 ms$ is driven by helical cooling. For example, the strong island growth in the O-point case is triggered by a strong cooling of the O-point for $t<1.1 ms$⁠, leading to an increase in the helical temperature perturbation (see Fig. 19, upper row). The perturbation decays primarily during the first crash. However, with growing island size, the current profile is modified in all cases, which may have a stabilizing or destabilizing effect for the continued island evolution. In the NPI-case, the current density gradient slightly inside the q = 2 surface remains negative, thus acting in a destabilizing manner throughout the first crash (see Fig. 20). In contrast, in the O-point case, a stabilizing, positive current gradient at q = 2 leads, in conjunction with the decay of the helical perturbation, to a decrease of the island width during the first crash. However, the gradient changes sign around $t=1.9 ms$⁠, subsequently leading to a destabilization.

In the X-point case, initially (⁠$t=0.6,…,1.0 ms$⁠) a cooling of the X-point reduces the preexisting temperature perturbation (see Fig. 19, bottom row), which causes the observed slight decrease in the island width. A distinct perturbation becomes visible again after $t=1 ms$⁠, which rotates clockwise in the poloidal plane and drives the island growth and position. Simultaneously, the current gradient turns negative, which further destabilizes the island.

For the case with $ϕO=90°$⁠, a consistent behavior is observed: the current gradient also flattens during the first crash and only becomes negative around $t=1.7 ms$⁠. The island continues to grow from that point on.

A few additional simulations were performed that are briefly summarized in the following to assess specific details of the dynamics (we do not show figures or a more extensive discussion to avoid overly extending the length of the manuscript): the three runs described here were restarted from $t=1 ms$ with the configuration of $η=η0(Te,axi/Te,0)−3/2$⁠, introduced above in Sec. III A 1. This allows the linear drive to be viewed in isolation and to determine, if the bifurcating current profile evolution really causes the different behavior of the island dynamics. Indeed, the general trends remain the same: the 2/1 islands in the X-point and NPI-case continue to grow, albeit at reduced rates. In the O-point case, the island width remains constant and only decreases slightly after the first crash. It follows that helical effects nevertheless play a significant role for the island width reduction in the O-point case. In further tests, the simulations were repeated with smaller parallel heat conductivity. While a small reduction (by a factor of 1.6) does not greatly interrupt the dynamics, we see distinct differences for a reduction by a factor of 16: overall, the effect of the existing islands becomes smaller and the thermal quench in the O-point case is only delayed by about $0.2 ms$ compared to the NPI-case. During the first crash, the island growth rate stagnates only transiently, but we no longer see a reduction in size.

In summary, it can be stated that depending on the injection phase, the sign of the current gradient around the rational surface changes shortly after the injection in a way that either stabilizes or destabilizes the 2/1 mode during the first crash. In all cases, the gradient finally becomes negative sometime after the first crash, which causes further mode destabilization and consequently leads to the TQ. The exact origin of the different current profile evolution observed between the cases will require future work.

In the following, we describe the three cases with $NSPI=3×1020, NSPI=1×1021$⁠, and $NSPI=6×1021$⁠. First, we discuss these cases acting on an unperturbed plasma in Sec. III B 1 and compare it with the case with $NSPI=6×1020$⁠, which was discussed in detail in Sec. III A. The important cases of O-point injection and X-point injection into the large preexisting island (⁠$wi=6 cm$⁠) are described in succession in Sec. III B 2. For this large island size, the effect of the preexisting structures is again most pronounced. We find here that the X-point injection into a large preexisting islands tends to lower the threshold in the amount of material injected to trigger a TQ. Consequently, the X-point case with $NSPI=3×1020$ exhibits a TQ in contrast to the corresponding NPI-case.

For $NSPI=3×1020$⁠, i.e., the smallest amount of injected material studied here, full ablation is achieved 1 ms after the shards have started to affect the plasma, leading to a relative increase in the total particle content by about 150% corresponding to an assimilation of around 93% of the injected material. The core thermal confinement is affected significantly only from $t=1.7 ms$ on, i.e., several hundred microseconds after full ablation. Subsequently, after $t=1.9 ms$⁠, the decay of $Ecore$ abates, too. The core particle content has finally increased by 65% at around $t=2.5 ms$⁠. Full plasma stochastization is not reached within the simulated time such that the plasma would either recover or exhibit only a delayed TQ.

The case with $NSPI=1×1021$ behaves qualitatively similarly to the case with $NSPI=6×1020$⁠, discussed in Sec. III A 1: once again we see a core crash in the evolution of $Ecore$⁠, followed by a stabilization and a full TQ including full stochastization at $t=1.9 ms$⁠. However, the time between both crashes however is reduced to $0.5 ms$ and the material is fully ablated at $t=3.3 ms$⁠, i.e., fairly after the second crash. It can be compared to case C in the study³ (called H-mode case in this paragraph as opposed to the L-mode case studied here). In terms of SPI parameters, the setup is very similar and differs only by the number of shards (30 in contrast to 10 used here). We also observe a TQ shortly after injection in the H-mode case and a qualitatively similar evolution of the core temperature, which crashes about one millisecond after injection and within three hundred microseconds. The core density stays virtually unaffected in the H-mode case in contrast to the L-mode case due to a reduced penetration depth in the H-mode case because of the higher temperatures and the larger number of shards.

For the largest amount of material (⁠$6×1021$ atoms), the rapid drop of $Ecore$ begins as soon as $t=0.9 ms$ and $0.3 ms$ later, the total particle content is six times larger than the initial value, while only 1/6 of the injected material has been ablated. During the global loss of thermal confinement after $t=1 ms$⁠, local temperatures around the shards fluctuate, which leads to varying ablation rates. Full stochastization is reached around $t=1.2 ms$⁠. It should be noted that in runs with a large amount of material (⁠$NSPI=6×1021$⁠), the toroidal resolution may not be sufficient anymore and this could artificially induce an early full TQ.

To summarize the findings from this subsection and Sec. III A 1, we state that

• for $NSPI=3×1020$⁠, we see a temporary degradation of the heat confinement and an increase in the core particle content. However, the excited MHD activity is not sufficient to cause a prompt TQ.

• For $NSPI=6×1020$ or $NSPI=1×1021$⁠, we note the first impact on the core particle number at around $t=1.2 ms$ followed by a first burst of MHD activity. The plasma then stabilizes and a full stochastization only occurs with a delay, after a further $1.3$ or $0.7 ms$⁠, respectively. For $NSPI=6×1020$⁠, this is at a time point at which the shards have already fully ablated.

• Even more material triggers an immediate TQ at an early stage of the injection, when the shards are still far outside the plasma core and only a fraction of the material is ablated.

First, we discuss the O-point injections based on Fig. 21: we see a distinct impact by an island for the $NSPI=3×1020$ case. Apart from the intrinsic losses, $Ecore$ remains constant up to $2.4 ms$ and only then does the core crash occur, which is $0.7 ms$ later when compared to the NPI-case. It is worth noting that this crash happens clearly after completed ablation. In the further course of the process, $Ecore$ relaxes to the same magnitude as in the regarding NPI-case after about $t=3 ms$⁠. A full TQ does not occur, and the eventual appearance of this plasma is similar to the NPI-case in terms of particle and energy content and magnetic topology. In both cases, flux surfaces inside q = 1.5 remain fully intact. This behavior is an indication that injection parameters that are not sufficient to trigger a TQ can still not trigger it even in the presence of a big preexisting island for the case of O-point injection. This is in agreement with the observations made in Sec. III A 2 that O-point injection is to be favored for the proposed mitigation scheme.

Also, for the case with $NSPI=1×1021$⁠, the preexisting island acts in a retarding manner on the whole dynamics after the first crash and the full TQ occurs only at $t=2.4 ms$⁠, i.e., $0.6 ms$ later than in the NPI-case. For $NSPI=6×1021$⁠, the simulation suggests that there is no impact of the preexisting island and still a rapid TQ at $t=1.2 ms$⁠. The perturbations excited by the material injection dominate the preexisting perturbations at a very early stage.

We close this subsection with a look at the X-point injections, which are shown in Fig. 22: for $NSPI=3×1020, Ecore$ drops as early as $t=1.4 ms$⁠, which is $0.3 ms$ earlier and more violent compared to the corresponding NPI-case. After $t=2.2 ms$⁠, the dynamics deviate strongly from the NPI-case. Instead of it stabilizing, $Etot$ begins to drop rapidly. Accordingly, the plasma reaches full stochastization at around $t=2.2 ms$ and a TQ sets in. This is an remarkable observation: it indicates that a preexisting island can indeed reduce the amount of material needed to triggering a TQ when the injection phase coincides with the island X-point. Despite this potentially very important observation, no further analysis is made to prevent this article becoming too long.

To investigate how much the previous findings depend on the selected equilibrium, simulations were performed with a modified q-profile. To do so, the toroidal field amplitude was adjusted, to set q₀ to 1.1. Consequently, q = 2 is shifted outwards to $ΨN=0.76$⁠. Importing exactly the same perturbations as for the previous runs results in significantly smaller 2/1 island widths. For example, the perturbation, which led to $wi=6 cm$ in the cases discussed in Secs. III A and III B, here only results in an initial island width of $wi=4.5 cm$⁠. Simulations without preexisting islands show that $NSPI=6×1020$ does not trigger a full TQ in this equilibrium. A partial core crash can be identified at around $t=1.8 ms$ to $t=2.0 ms$⁠, where a 1/1 mode plays a key role. However, after $t=3.9 ms$ the n = 1 magnetic energy begins to decrease before reaching the critical amplitude required to provoke a full TQ.

The NPI-case, the X-point case, and the O-point case for $NSPI=6×1020$ and $wi=4.5 cm$ are shown in Fig. 23. Despite the greatly altered TQ dynamics and significantly smaller initial island widths, qualitatively similar mechanisms can still be observed. In the X-point case, the island width initially decreases between $t=0.5 ms$ and $t=0.8 ms$ and then starts to increase and continues to grow during the first crash. After the first crash (⁠$t=1.8 ms$⁠), the dynamics are similar to the NPI-case. In the O-point case, the island begins to grow initially at $t=0.5 ms$⁠, stagnates at around $t=1.1 ms$⁠, and then slightly decreases in the following $0.2 ms$⁠. Note that this happens clearly before the first crash. The growth after the crash is now clearly delayed by $0.3 ms$⁠, and the saturated island width is only reached after $t=4.2 ms$⁠. This saturated island width is of similar amplitude to that in the NPI-case. These observations are indeed in agreement with the cases investigated previously: when we compare the X-point and O-point cases, we can identify a first sequence, where shards are closer to the q = 2 surface in which $w2/1$ is larger in the O-point case. After a turning point (here around $t=1.5 ms$⁠), the island in the X-point case is now always the larger one. Comparing with the NPI-case confirms, that the injection into the O-point still has a delaying effect on the long-term evolution of the island. Nevertheless, injecting into the X-point does eventually lead to a larger saturated island width. This is in agreement with the observations made for the other q-profile and X-point injection with $NSPI=3×1020$⁠.

As these are only first preliminary findings for the modified q-profile, more studies are required. This involves further scanning of $NSPI$ to find a value for a real TQ triggering, as well as importing of larger perturbations to study the effect of larger (e.g., $wi=6 cm$⁠) preexisting islands.

To analyze the applicability of massive deuterium injection for fast plasma dilution in the presence of preexisting tearing modes and, in particular, the 2/1 NTM, a series of simulations were performed. To do so, an ASDEX Upgrade L-mode-like plasma with an injecting nozzle on the outer midplane was considered. The simulations included realistic temperature dependent plasma resistivity and parallel heat conductivity as well as an evolution of the bootstrap current fraction. Background impurities and Ohmic heating were ignored, and the initial central safety factor was well above unity. Ten shards were considered, representing a broken pellet, from which the most effective penetration properties can be expected. Scans over the amount of material injected, the amplitude of the preexisting perturbations, and the relative toroidal phase between the nozzle and the preexisting island were carried out.

Scanning over the amount of material injected into a plasma without preexisting islands showed that for $NSPI=6×1020$⁠, the shards trigger a partial crash around $1.1 ms$ after the injection (at $t=0.1 ms$⁠). From this crash on, the core is fueled steadily. At $t=1.9 ms$⁠, shards are close to the magnetic axis, where they become fully ablated. Only at $t=2.5 ms$⁠, does the TQ occur. This case is particularly interesting for this study, as it can be understood as a limiting case between the two regimes of no TQ and an immediate TQ, prior to full ablation. Indeed, a significant effect on the TQ timing could be identified for larger preexisting islands. If injecting with $NSPI=6×1020$ into the O-point of an island of initial width $wi=6 cm$⁠, the TQ is delayed by about $1 ms$⁠. This behavior, which appear to be favorable for core dilution, is caused by shrinking of the island width after the first crash up to $1.9 ms$⁠. However, injecting directly into the X-point results in the opposite reaction. Rapid island growth during the first crash causes a full TQ as early as $t=1.8 ms$⁠. By comparing the island size evolution of both cases, two sequences are identifiable: up to $t=1.0 ms$⁠, the island size is larger when injecting into the O-point. However, injecting around the X-point has a temporary reducing effect onto the island width. By contrast, the subsequent evolution of the current density profile in both cases differs significantly, leading to a strong destabilization for the case of X-point injection and a temporary stabilization in the O-point case. As a result, the island in the X-point case becomes larger after $t=1 ms$⁠, while, in the O-point case, the island growth reestablishes only at $t=1.9 ms$⁠. Helical cooling inside the 2/1 island has long been identified as an important drive for tearing modes and a key player in disruptions triggered by massive material injection; however, our results suggest that other mechanisms are also involved. Indeed, if helical cooling was the only mechanism involved, we would have expected to observe a faster TQ triggering with the O-point injection, whereas we observe the opposite. The additional mechanisms may be related to modifications to the pressure distribution caused by the injection.

In respect to the main questions, these series of simulations show that preexisting islands do not seem to preclude the option for dilution by massive deuterium injection as part of a disruption mitigation strategy. While larger islands have a significant impact on TQ timing, the effect of smaller islands is approaching negligible. In all cases, dilution is scarcely affected by the preexisting prior to the TQ. Furthermore, these simulations indicate that the TQ triggering threshold in the amount of injected deuterium is only decreased if it is injected around the X-point. This was observed for the case with $NSPI=3×1020$⁠, which does not trigger a TQ if applied to a healthy plasma, but where a TQ is seen if injecting into the X-point of the largest considered preexisting island. These simulations further suggest that injecting around the X-point of a big preexisting island should generally be avoided, as this may have a premature effect on the TQ timing.

These first findings need to be validated by taking into account other physical mechanisms as well as ITER scenarios. Background impurities can have an amplifying effect on island growth if they further cool down the O-point by radiation, thereby enhancing helical perturbations. Their role for the considered strategy will be assessed in future studies for an ASDEX Upgrade. In addition, the exact processes in a case with q₀ close to unity need to be investigated in more detail. However, the preliminary results show qualitative agreement between the equilibria with different q-profiles. In this study, only one shape of a (very large) neutral cloud was considered. The role of size and shape of the cloud requires further investigation. In particular, a possible relationship between the cloud size and the minimum width of an island required to cause a measurable impact could be investigated. In further steps, this analysis will be extended to ITER. In larger scale machines, larger time scales of TQ triggering would be of benefit for the considered strategy, while the material will also take more time to penetrate and dilute the plasma thoroughly. Multiple injections like those planned for ITER will also be a focus of further studies. As the shard clouds will typically not arrive at the plasma exactly at the same time, it is expected that the interaction of a “leading shard cloud” with the preexisting island will dominate and the overall mechanisms of a multiple injection scheme will remain very similar to what has been studied here.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training Program Nos. 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Some of the work was performed using the Marconi-Fusion supercomputer.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.