Numerical study of runaway current impact on sawtooth oscillations in tokamaks

1.   Introduction

The presence of high-energy runaway electrons (REs) poses a significant threat to the successful operation of large tokamaks, especially as the plasma current ($I_{\text{p}}$) increases. The likelihood of generating high-current RE beams escalates with rising $I_{\text{p}}$ due to the exponential sensitivity of the RE avalanche multiplication process to this parameter [1]. This risk is particularly pronounced in fusion reactors like ITER during disruptions, often resulting in the activation and melting of plasma-facing components due to the loss of a RE beam with energies up to tens of MeV and currents up to about 10 MA [2–4]. Therefore, it is essential for reactor-scale tokamaks to either prevent RE generation or enable the rapid dissipation of increasing RE currents to safely achieve their intended $I\mathrm{_p}$.

Recent efforts to understand and mitigate RE effects for ITER on existing tokamaks involve interpreting and modeling experimental measurements with innovative diagnostics, as well as comparing simulation results with observed phenomena [1, 5–18]. Various strategies have been employed, including injecting materials into the plasma [7, 19–22] to enhance density and collisions, thereby increasing the friction force, and introducing 3D magnetic field perturbations externally [23–25] to diffuse REs out of the plasma before they reach higher energies and undergo avalanche multiplication. In addition, the degradation of RE confinement by MHD instabilities is well-known [2]. Thus, another passive approach involves utilizing macroscopic MHD instabilities of various kinds to dissipate the RE beam [17, 18, 26–30]. For instance, it has been proposed that sudden impurity injection during the current quench (CQ) could destabilize MHD due to current profile shrinking, thereby enhancing RE loss to the wall sufficiently to counteract the CQ avalanche gain.

Previous studies have explored the interplay between REs and resistive MHD instabilities in post-disruption plasmas [14, 15, 17, 18, 31]. Experimental evidence from DIII-D has shown sawtooth-like relaxation of the RE current profile but without triggering RE loss in a post-thermal quench (TQ) RE beam [13]. Moreover, modeling of plasma equilibrium evolution during a disruption for ITER has revealed cases where the on-axis safety factor $q_0$ drops below unity [32]. These findings, combined with the low temperature of the post-TQ plasma for the bulk species, suggest that a resistive internal kink could become unstable, potentially producing 3D perturbations that may de-confine REs in the beam [18]. However, a simulation study employing M3D found that sawtooth oscillations in the runaway plasma were suppressed due to the flattening of the current profile after a single occurrence of the resistive internal kink mode [33]. In addition, this study did not examine the subsequent loss of REs. Given the significance of this issue, further in-depth investigation into the dynamics of sawtooth oscillations in the presence of REs and the resulting RE loss is essential, as these MHD instabilities hold promise in mitigating RE damage in a tokamak reactor. Fortunately, the integration of a fluid model for the RE current into the fully 3D MHD codes M3D [34], M3D-C$^1$ [35, 36], JOREK [37, 38] and NIMROD [39], has enabled the study of MHD instabilities, including resistive internal kink, tearing and hose modes [14–18, 33, 40], which thus serves as the primary motivation for this work.

This work aims to examine the internal kink instability effected by REs and its impact on RE loss using an implicit extended-MHD code M3D-C¹. The results indicate that RE can destabilize the resistive internal kink and induce non-zero real frequencies of the modes, which are consistent with those obtained by MASR-F [17] and M3D-C¹ [15, 17]. In the nonlinear phase, with the assumption that the runaway current entirely replaces the plasmas current, it is noted that the sawtooth oscillations are significantly suppressed. More specifically, the first sawtooth cycle occurs earlier compared to the case without runaway current. After the first sawtooth collapse, the safety factor $q$ slightly exceeds unity within the $q=1$ surface. Meanwhile, both plasma current density and runaway current density flatten within this surface, albeit displaying fine structures. Subsequently, the high-$n$ and high-$m$ modes, emerging from the nonlinear couplings of low-$n$ modes, become much more unstable in the presence of REs, with $n$ and $m$ being the toroidal and poloidal mode numbers, respectively. The growing fluctuations of these high-$n$ and $m$ modes impact the saturation level of the resistive internal kink mode, particularly resulting in several small drops in core electron temperature following the first sawtooth oscillation. The study also shows that despite considerable 3D perturbations linked to internal kink and high-$n$ and high-$m$ modes, REs stay well-confined in the plasma core region (within the $q=1$ surface). In addition, the final runaway current density profile shows a heightened steepening near the $q=1$ surface and a distinctive corrugated inhomogeneity within it when compared to the initial profile, attributed to the shrinking of the current channel and the improved confinement of REs within $(1,1)$ magnetic islands, along with the impact of the increasingly growing $n=0$ component fluctuation.

The paper is organized as follows. Section 2 provides a brief description of the fluid RE model implemented in the M3D-C$^1$ code. Section 3 discusses the simulation results of a resistive internal kink instability influenced by REs, focusing on the mode growth rate, real frequency, saturation level, and subsequent RE loss. In addition, a detailed comparison of our findings with those presented in [33] is presented, highlighting the significance of high-$n$ and high-$m$ mode dynamics on sawtooth oscillations—dynamics that were overlooked in their study. Finally, section 4 concludes with a summary of the research.

2.   Simulation model and setting

We employ the implicit MHD code M3D-C¹ [35, 36], which utilizes high-order C¹ continuous finite elements in three dimensions, where C¹ denotes continuity in first derivatives. The code is capable of simulating various physical models, including ideal-MHD, resistive-MHD, and two-fluid models. In this work, we solve the following set of MHD equations with the recently implemented fluid model of REs [14, 15, 18].

Here, equation (1) is the MHD momentum equation, where $n$ is the plasma number density, $M_\text{i}$ is the ion mass, ${{{\boldsymbol{V}}}}$ is MHD velocity, ${{{\boldsymbol{E}}}}$ is the electric field, $e$ is the elementary charge, ${{{\boldsymbol{J}}}}$ and ${{{\boldsymbol{J}}}}_\text{RE}$ are, respectively, the total current and the runaway current, ${{{\boldsymbol{B}}}}$ is the magnetic field, and $p$ is the total plasma pressure. The first term on the right-hand-side accounts for the positive charge of plasma, excluding RE. Equation (2) is the continuity equation for plasma density, while equation (3) is the continuity equation for RE density $n_\text{RE}$, with ${{{\boldsymbol{b}}}}={{{\boldsymbol{B}}}}/B$, $C_\text{RE}$ denoting the RE convection velocity, and $S _ \text{RE}$ as a source term for RE generation. Note that here, only the parallel velocity $C_\text{RE}$ and ${{{\boldsymbol{E}}}}\times{{{\boldsymbol{B}}}}$ drift are considered, neglecting magnetic gradient and curvature drifts for the sake of simplicity. However, these magnetic drifts of the REs can be important in the context of equilibrium and MHD behavior. The extent to which neglecting these magnetic drifts can affect the solutions is not yet fully understood and will be investigated in future work. Typically, the value of $C_\text{RE}$ is chosen slightly below the speed of light $c$ for computational efficiency, as REs are mainly relativistic particles with very small pitch angles [18]. Furthermore, it should be emphasized that the generation source of REs is not considered in this study ($S _ \text{RE}=0$). Given that the time scale of MHD instabilities is significantly shorter than that of RE generation, MHD instabilities can be examined with a prescribed initial runaway current profile without a runaway source.

The RE current ${{{\boldsymbol{J}}}}_\text{RE}$ is derived from $n_\text{RE}$ in equation (4). Upon substitution of ${{{\boldsymbol{J}}}}_\text{RE}$ into equation (1), the term of perpendicular electric field is cancelled. Nevertheless, the current implementation omits the remaining parallel component due to the significantly lower RE density compared to thermal electrons and the relatively weaker ${{{\boldsymbol{E}}}}$ parallel to ${{{\boldsymbol{B}}}}$ compared to the perpendicular component [15, 18, 31, 33]. Given the negligible collisional friction of REs compared to thermal electrons, the resistive term in the generalized Ohm’s law is simplified, as shown in equation (5), where $\eta$ represents the plasma resistivity. Equation (6), Faraday’s law, is utilized to update the vector potential ${{{\boldsymbol{A}}}}$ to maintain a divergence-free magnetic field ${{{\boldsymbol{B}}}}={\boldsymbol\nabla}\times{{{\boldsymbol{A}}}}$. Equation (7) is the pressure equation, where $\gamma$ is the ratio of specific heats and ${{{\boldsymbol{q}}}}$ is the heat flux. The last term in this equation accounts for Ohmic heating. The heat flux ${{{\boldsymbol{q}}}}$ can be calculated with equation (8), where $\kappa_\parallel$ and $\kappa_\perp$ denote parallel and perpendicular heat conductivities. Note that the equilibrium field $n^0$, ${{{\boldsymbol{J}}}}^0$, and $T^0$ subtracted in the dissipation terms act as effective current and heat sources, which are necessary in simulating sawtooth oscillations.

The equations of M3D-C¹ are formulated in cylindrical $(R,\phi,Z)$ coordinates. The split time advance scheme in M3D-C¹ involves computing the MHD equations (8)−(1) for each MHD timestep using the $\theta-$implicit method for time-integration. The momentum equation is advanced semi-implicitly, incorporating a parabolic term derived from second-order time derivatives to ensure numerical stability [35]. Subsequently, plasma and RE densities are calculated using magnetic and electric fields from the previous timestep. The magnetic field and the pressure equations are advanced last, using velocities and RE density from both the current and previous timesteps.

In M3D-C¹, the system of units for all quantities presented hereafter is derived from a characteristic length, density, and magnetic field. These normalizations are based on a normal length $L_0=1$ m, a normal particle density $n_0=1\times 10^{20}$ ${\rm m}^{-3}$ for a hydrogen plasma, and a normal strength of the magnetic field $B_{0}=1$ T. This establishes the normalization for velocity as $\upsilon_\text{A0}=B_{0}/(\mu_0n_0M_\text{i})^{1/2}$ the Alfvén speed, the time to $\tau_\text{A0}=L_0/\upsilon_\text{A0}$, pressure as $p_0= B_0^2/\mu_0$, electric field as $E_0=\upsilon_0B_0/c$, current as $I_0=B_0cL_0/\mu_0$, and current density as $J_0=cB_0/\mu_0L_0$. In addition, resistivity $\eta$, mass diffusivity $D$, perpendicular and parallel thermal conductivities $\kappa_\perp$ and $\kappa_\parallel$, and viscosity $\nu$ are normalized as follows: ${\hat \eta}\equiv\eta/(\mu_0L_0^2/\tau_\text{A0}c^2)$, ${\hat D}\equiv D/(L_0^2/\tau_\text{A0})$, ${\hat \kappa}_\perp \equiv \kappa_\perp/(n_0L_0^2/\tau_\text{A0})$, ${\hat \kappa}_\parallel \equiv \kappa_\parallel/(n_0L_0^2/\tau_\text{A0})$, and ${\hat \nu}\equiv \nu/(B_0^2\tau_\text{A0}/\mu_0)$, respectively.

We examine a toroidal limiter plasma with a major radius of $R_0=1.65$ m and a minor radius of $a=0.4$ m, featuring an approximately circular poloidal cross-section. The on-axis magnetic field strength is $B_\text{t}=1.28$ T. These parameters are similar to the HL-2A tokamak, but our analysis emphasizes the physical mechanism rather than specific devices. The initial toroidal current density $J_\phi$ and safety factor $q$ profiles are shown in figure 1(a). The safety factor profile $q$ is from the EFIT reconstruction [41]. The current density profile self-consistently satisfies the Grad-Shafranov equation, with an on-axis safety factor of $q_0=0.8$ for an unstable internal kink mode. The $q=1$ surface is at $\psi_\text{N}=0.20$ with $\psi_\text{N}$ the normalized poloidal flux defined as $\psi_\text{N}=(\psi-\psi_0)/(\psi_{\text b}-\psi_0)$, where $\psi_0$ and $\psi_\text{b}$ are the poloidal flux at the magnetic axis and at the plasma boundary, respectively. Figure 1(b) shows the finite element mesh used in the linear simulations, which is refined near the $q=1$ flux surface. For 3D nonlinear cases, an unrefined mesh with 8 toroidal planes and over 5000 elements per plane is used. The parameters in normalized units of M3D-C$^1$ adopted here are ${\hat\eta}\simeq 1.5\times 10^{-6}$, $C_\text{RE}/\upsilon_\text{A0}=50$, timestep $\Delta t/\tau_\text{A0}=5.0$, ${\hat\nu}= 2.0\times 10^{-6}$, ${\hat D}=1.0\times 10^{-6}$, ${\hat\kappa}_\perp=2.0\times 10^{-6}$, and ${\hat\kappa}_\parallel=10^6{\hat\kappa}_\perp$.

Figure  1.  (a) The initial profiles of safety fator $q$ and toroidal current density $J_\phi$ as functions of the normalized poloidal flux $\psi_\text{N}$; (b) finite element mesh used in M3D-C$^{1}$ simulation. Shown are the finite element mesh (bule), the $q=1$ surface (red) and the last closed flux surface (purple). The mesh has approximately 1.1$\times$ $10^4$ elements.

DownLoad: Full-Size Img PowerPoint

3.   Simulation results and discussions

In this section, we first examine the effects of REs on the linear properties of the resistive internal kink mode. Figure 2 shows the dependence of (a) growth rate ($\gamma$) and (b) real frequency ($\omega$) for the $(m,n)=(1,1)$ resistive internal kink mode, comparing cases with and without runaway current. It is important to emphasize that, for the case with RE current discussed in this study, we assume that all equilibrium current is carried by REs; that is, ${{{\boldsymbol{J}}}}^0={{{\boldsymbol{J}}}}^0_\text{RE}$, which is characteristic of the post-disruption RE plateau. In the absence of runaway current, depicted in blue in figure 2(a), the growth rate deviates from the expected $1/3$ power law for high $\eta$ values, consistent with findings reported in [42]. Upon incorporating the RE effect, the growth rate shown in red exhibits an improved recovery of the $\eta^{1/3}$ scaling at higher $\eta$ values, in agreement with results from the MARS-F code as noted in [17]. Furthermore, figure 2(b) highlights that the presence of REs introduces a finite mode frequency, which increases with increasing plasma resistivity. This observation is consistent with theoretical predictions and simulations detailed in earlier studies [14, 15, 17].

Figure  2.  Scalings of the (a) growth rate $\gamma$ and (b) real frequency $\omega$ of $(1,1)$ resistive kink mode versus resistivity in linear simulations, for cases with ${{{\boldsymbol{J}}}}_\text{RE}=0$ (blue curve with circle markers) and with ${{{\boldsymbol{J}}}}^0_\text{RE}={{{\boldsymbol{J}}}}^0$ (red curve with square markers).

DownLoad: Full-Size Img PowerPoint

We now investigate the impact of REs on the nonlinear evolution of the mode, starting with the results of the sawtooth cycles in the absence of RE effects. Figure 3 shows the time evolutions of (a) total kinetic energy, (b) core electron temperature $T_\text{e;core}$ measured at $(R,\phi,Z)=(1.658,0,0)$ with $T_\text{e}=T/2$, and (c) kinetic energy for different toroidal mode numbers. Typical sawtooth cycles are observed, with the first cycle influenced by the initial conditions. Over subsequent cycles, the system establishes a recurring pattern. In figure 3(c), it is clear that the most unstable mode is dominated by the $n=1$ component. Figure 3(b) shows that $T_\text{e;core}$ periodically peaks until a rapid onset of instability occurs at the $q=1$ surface, leading to magnetic reconnection and a sudden temperature flattening near the center, as highlighted by the blue and red curves in figure 4. Following a crash, the $q$ profile flattens within the core region. In this process, the subtracted equilibrium fields act as effective sources to restore the initial equilibrium, enabling the $q$ profile to drop below unity at the magnetic axis (as shown by the red curve in the zoomed-in figure of figure 4(b)), triggering the subsequent sawtooth event. In addition, the Poincaré plots of magnetic field lines during one sawtooth cycle, marked by the vertical dashed lines in figure 3(b), are shown in figure 5, which confirms the periodic emergence and disappearance of the $(1,1)$ magnetic island.

Figure  3.  The time evolutions of (a) total kinetic energy, (b) the core electron temperature $T_\text{e;core}$ and (c) kinetic energy of the $n=0$ –4 components for sawtooth in the absence of runaway electrons.

DownLoad: Full-Size Img PowerPoint

Figure  4.  (a) The changes in electron temperature profiles for the selected time slices; (b) the dependence of $q$ on $\psi_\text{N}$ before, during and after the sawtooth collapse, where the black, blue and red lines correspond to the instants indicated by the vertical dashed lines in figure 3(b).

DownLoad: Full-Size Img PowerPoint

Figure  5.  Poincaré plots of the magnetic field lines during one sawtooth cycle shown in figure 3 in the absence of REs. (a) $t=4100\tau_\text{A0}$, (b) $t=4525\tau_\text{A0}$, (c) $t=4750\tau_\text{A0}$, and (d) $t=5000\tau_\text{A0}$.

DownLoad: Full-Size Img PowerPoint

Next, we examine the influence of REs on sawtooth oscillations. Figure 6 shows the time evolutions of (a) total kinetic energy, (b) $T_\text{e;core}$, and (c) kinetic energy for the $n=0$ –4 components during sawtooth oscillations in the presence of REs. When the plasma current is entirely carried by runaway current, the sawtooth oscillations are significantly suppressed due to the absence of resistive diffusion in the runaway current, thereby preventing its return to the original profile. It is noteworthy that the value of $q$ at the magnetic axis ($q_0$) does not consistently remain above unity after the first crash around 1000$\tau_\text{A0}$. Subsequently, three perturbations in temperature occur between 2350$\tau_\text{A0}$ and 4600$\tau_\text{A0}$, highlighted by distinct color-coded phases in the zoomed-in section of figure 6(b), which will be discussed in detail later. Throughout the entire nonlinear phase, the $n=1$ mode remains dominant, as shown in figure 6(c). After the first crash (approximately $t=1800\tau_\text{A0}$), besides the $n=1$ mode, the high-$n$ modes, predominantly the $n=2$ mode depicted in the zoomed-in section of figure 6(c), start to grow, potentially affecting the saturation level of the $(m,n)=(1,1)$ mode and the confinement of REs.

Figure  6.  The time evolutions of (a) total kinetic energy, (b) the core electron temperature $T_\text{e;core}$ and (c) kinetic energy of the $n=0$ –4 components for sawtooth in the presence of runaway electrons.

DownLoad: Full-Size Img PowerPoint

To delve into the details of the first sawtooth cycle, we analyze three different time slices: before, during, and after the crash, as indicated by the vertical dashed lines in figure 6(a). Figure 7 shows the profiles of toroidal current density ($J_{\phi}$), runaway current density ($J_\text{RE}$), and the safety factor ($q$). After the first crash, it becomes apparent that both $J_{\phi}$ and $J_\text{RE}$ along with $q$, are significantly flattened within the $q=1$ surface compared to their pre-crash profiles. Furthermore, the value of $q$ in the core slightly exceeds unity. These results are consistent with those reported in reference [33], although a notable distinction is the presence of fine structure in the flattened profiles, which will be discussed later. Note that, during the crash, there is a peak in $J_{\phi}$, which leads to instability and subsequent reconnection, ultimately resulting in the broadening of the current profile.

Figure  7.  Profiles of (a) toroidal current density, (b) runaway current density, and (c) safety factor ($q$) at three different time slices: $t=750\tau_\text{A0}$ (before the crash, shown by the black curve), $t=1125\tau_\text{A0}$ (during the crash, shown by the blue curve), and $t=1500\tau_\text{A0}$ (after the crash, shown by the red curve). This case involves the presence of REs.

DownLoad: Full-Size Img PowerPoint

Following the first crash, three subsequent smaller $T_\text{e;core}$ drop. Figures 8−10 present the selected Poincaré plots at various time slices for each phase, illustrating the evolution of the magnetic field. In the phase marked as ①, the region near the magnetic axis experiences $(1,1)$ resistive internal kink instabilities due to $q_0<1$, resulting in typical sawtooth periodic behavior. While the magnetic surfaces inside remain relatively intact, the growing magnetic islands of high-$m$ modes, such as $(2,1)$ and $(3,1)$, lead to increasing stochasticity in the magnetic surfaces outside, as illustrated in figure 8. In the phase marked as ②, the area near the magnetic axis continues to display typical sawtooth periodic behavior. However, the magnetic surfaces outside exhibit complete stochastic characteristics during the oscillation period, as shown in figure 9. In the phase marked as ③, magnetic reconnection occurs at $q_\text{min}$ with $\psi_\text{N}=0.2175$ in a region of weak reverse magnetic shear, while the magnetic surfaces inside the $q=1$ surface remain intact. Compared to phase ②, the magnetic surfaces outside show signs of restoration, as demonstrated in figure 10. By comparing figures 9 and 5, we note that in the presence of REs, the width of the (2,1) and (3,1) magnetic islands increases compared to the case without REs. This increase suggests that these high-$m$ modes could become unstable due to the steepening of the current density profile at their rational surfaces. The redistribution of current density is attributed to the effect of REs on the instability of the (1,1) resistive kink mode, which in turn affects the dynamics of the nonlinear mode couplings. This is further illustrated in figures 13(b), where the profiles of runaway current density display fine structure in the radial direction.

Figure  8.  Poincaré plots of the magnetic field lines during the period marked as ① in figure 6(b), corresponding to the case involving the presence of REs. (a) $t=2350\tau_\text{A0}$, (b) $t=2550\tau_\text{A0}$, and (c) $t=2750\tau_\text{A0}$.

DownLoad: Full-Size Img PowerPoint

Figure  10.  Poincaré plots of the magnetic field lines during the period marked as ③ in figure 6(b), corresponding to the case involving the presence of REs. (a) $t=4175\tau_\text{A0}$, (b) $t=4375\tau_\text{A0}$, and (c) $t=4575\tau_\text{A0}$.

DownLoad: Full-Size Img PowerPoint

Figure  9.  Poincaré plots of the magnetic field lines during the period marked as ② in figure 6(b), corresponding to the case involving the presence of REs. (a) $t=3000\tau_\text{A0}$, (b) $t=3150\tau_\text{A0}$, and (c) $t=3325\tau_\text{A0}$.

DownLoad: Full-Size Img PowerPoint

Figure  13.  The profiles of (a) the safety factor $q$ and (b) the runaway current density $J_\text{RE}$ are compared at $t=0$ (initial state shown by black curve), $t=2375\tau_\text{A0}$ (during crash shown by green curve) and $t=5000\tau_\text{A0}$ (after crash shown by purple curve, with the corresponding $q\leq 1$ region highlighted in yellow) for second sawtooth collapse with runaway current.

DownLoad: Full-Size Img PowerPoint

It is important to note that in the study by Cai et al [33], the influence of high-$n$ modes on saturation level of the (1,1) mode was overlooked. This could lead to the complete suppression of sawtooth oscillations after a single cycle. To validate this, we replicated the results of [33] by reducing the number of toroidal planes $n_\text{p}$ from the original 8 to 4. These toroidal planes represent the number of finite elements loaded toroidally. The simulation results presented in figure 11 are consistent with the findings of [33]. By combining the convergence analysis illustrated in figure A1(a) in the Appendix A with the results in figure 11, it is evident that after $t=1800\tau_\text{A0}$, both the $n=1$ and $n=2$ modes exhibit significant growth. A Fourier spectral decomposition of perturbed stream-function, $\delta U_{mn}$ [43], shown in figure A3 in Appendix B, identifies the contributing poloidal harmonics and indicates that the growth of the $n=1$ mode is dominated by the $(m,n)=(2,1)$ mode. For the $n=2$ mode, the dominated poloidal harmonics correspond to the $(m,n)=(4,2)$ mode. These findings indicate that the high-$n$ and high-$m$ modes significantly affect on the saturation level of the $(1,1)$ mode. In addition, the Poincaré plots in figure 12 for selected time slices indicated by the vertical lines in figure 11(b) exhibit relatively weak stochasticity compared to those in figures 8−10. This suggests that the high-$n$ modes are much more unstable in the presence of REs, playing a crucial role in the formation of stochastic magnetic fields. We remark here that the growing high-$n$ and high-$m$ modes could be driven by a steepening of the current density profile at their rational surfaces, as well as by the nonlinear couplings of low-$n$ modes.

Figure  11.  The time evolutions of total kinetic energy for sawtooth (a) in the absence of runaway electrons and (b) in the presence of runaway electrons in the case of $n_\text{p}=4$.

DownLoad: Full-Size Img PowerPoint

  A1.  The time evolutions of kinetic energy of the $n=0-3$ components for sawtooth in the presence of REs with different toroidal plane numbers.

DownLoad: Full-Size Img PowerPoint

Figure  12.  The evolution of the magnetic Poincaré plots for the selected time slices shown in figure 11(b), corresponding to (a) $t=1150\tau_\text{A0}$, (b) $t=2175\tau_\text{A0}$, (c) $t=2800\tau_\text{A0}$, and (d) $t=5000\tau_\text{A0}$. This case involves the presence of REs with the number of toroidal planes $n_{\rm p}=4$.

DownLoad: Full-Size Img PowerPoint

Then we specifically analyzed two time slices after the first crash, marked by the green and purple vertical lines in figure 6(a). A comparison of the safety factor ($q$) and runaway current density profiles ($J_{\rm RE}$) at $t=2375\tau_\text{A0}$ and $t=5000\tau_\text{A0}$ after the first sawtooth collapse is shown in figure 13. Notably, at $t=2375\tau_\text{A0}$, $q_0$ marginally drops below 1 but exceeds 1 again at $t=5000\tau_\text{A0}$, as shown in zoomed-in view of figure 13(a). In figure 13(b), the profile of $J_{\rm RE}$ after the crash (purple curve) shows a significantly steeper gradient near the $q=1$ surface compared to its initial profile. This trend is further illustrated by the surface plots of $J_{\rm RE}$ in the ($R, Z$) plane in figure 14, indicating a transition in the distribution of $J_{\rm RE}$ from a Gaussian-like shape to one characterized by increased steepness and corrugated inhomogeneities. This enhanced steepness is attributed to the shrinking of the current channel and improved confinement of REs within the $(1,1)$ magnetic islands, as depicted in figures 10 and 14. Moreover, the observed corrugated profiles are influenced by the $n=0$ component, with the corrugation becoming more pronounced as this component significantly grows, evidenced by the black curves in figure A1. These corrugated inhomogeneities may link to observations from DIII-D, captured by a fast visible camera and presented in figure 8 of reference [13], which reported that the flashes of visible radiation localized within the core of the RE beam evolved from solid circles into bright rings.

Figure  14.  Surface plots of the runaway current density ($J_{\rm RE}$) in the ($R$, $Z$) plane at four different time slices: (a) $t=750\tau_\text{A0}$ (indicated by the black vertical line in figure 6(a)), (b) $t=1125\tau_\text{A0}$ (indicated by the blue vertical line in figure 6(a)), (c) $t=2375\tau_\text{A0}$ (indicated by the green vertical line in figure 6(a)), and (d) $t=5000\tau_\text{A0}$ (indicated by the red vertical line in figure 6(a)). This case involves the presence of REs.

DownLoad: Full-Size Img PowerPoint

Finally, we examine the impact of sawtooth oscillations on RE loss. The Poincaré plots in figures 9 and 10 illustrate the emergence of $(2,1)$ and $(3,1)$ magnetic islands after the first sawtooth collapse cycle. Simulation results indicate that the growing high-$n$ and high-$m$ modes disrupt the magnetic surfaces, leading to stochasticity in the magnetic field lines in the outer region and consequently resulting in the loss of REs outside the $q=1$ surface. However, this effect appears to have minimal impact on the majority of REs confined within the $q=1$ surface, as demonstrated by the time evolution of the runaway current shown in figure 15. The overall fraction of runaway current lost during this period is approximately 3%. These findings are consistent with the experimental observations in DIII-D [13], where sawtooth-like relaxation of the RE current profile did not result in RE loss post-thermal quench. In figure 15, the intervals marked by the blue and red vertical lines encompass phases ① and ②, as shown in figure 6(b). It is noteworthy that, in addition to the rapid decline in runaway current during these phases, the subsequent decrease can be attributed to two aspects. First, the residual REs outside the $q=1$ surface, indicated beyond the yellow region in figure 13(b), continue to be lost due to the stochastic evolution of the magnetic surfaces. The second one may involve the weak reconnection of the $(1,1)$ mode at the $q=1$ rational surface, as shown in figure 10, which enhances the transport of central REs towards the tokamak edge through the stochasticity of magnetic field lines beyond the $q=1$ surface.

Figure  15.  The time evolution of the runaway current.

DownLoad: Full-Size Img PowerPoint

4.   Conclusion and outlook

This study have investigated the interplay between relativistic runaway electrons (REs) and the $n=1$ resistive internal kink instabilities in toroidal plasmas, particularly in the post-thermal quench, by employing the advanced capabilities of the 3D MHD code M3D-C$^{1}$. Our results show that the presence of REs significantly alters the behavior of the resistive internal kink mode, indicating their ability to destabilize this mode and induce non-zero real frequencies, consistent with theoretical predictions and earlier simulations. Notably, the transition to the nonlinear phase reveals a clear reduction in sawtooth oscillations, with the first sawtooth event occuring earlier when runaway currents are present. This early occurrence suggests significant changes in how the plasma responds to disruptions.

The flattening of both plasma and runaway current density profiles within the $q=1$ surface after the sawtooth collapse indicates important changes in confinement dynamics. We observe that these changes are characterized by fine structural features influenced by the $n=0$ component, and small temperature drops in the core are linked to the activation of high-$n$ and high-$m$ modes due to REs, which play a key role in the stochastic behavior of magnetic surfaces, thereby affecting the saturation level of the $(1,1)$ resistive internal kink mode.

Importantly, despite the considerable 3D perturbations associated with the resistive internal kink dynamics and the destabilization of high-$n$ and high-$m$ modes, our findings indicate that REs remain well-confined in the plasma core, particularly in the $q=1$ surface. The resulting runaway current density profile shows a notable steepening near this surface, along with distinct corrugated inhomogeneities, which can be attributed to the shrinking of the current channel and improved confinement within the $(1,1)$ magnetic islands.

We have also validated our findings through convergence testing, confirming the reliability of our simulation parameters and providing deeper insights into the dynamics of the $n=1$ mode and the behavior of REs in the post-thermal quenching phase. Looking ahead, future research will focus on studying the thermal and current quench phases at the onset of disruptions, particularly regarding the generation of REs in a self-consistent manner. By enhancing our understanding of plasma behavior during these disruption events, we aim to inform the design of more effective strategies for managing and mitigating similar phenomena in next-generation fusion reactors.

Appendix A.   Convergency tests

In this appendix, convergence tests are conducted on simulation parameters to ensure the accurate representation of the $n=1$ mode. Parameters such as the time step $\Delta t$, the RE convection velocity $C_\text{RE}$, and the number of toroidal planes $n_\text{p}$ are validated. Figure A1 shows the time evolutions of kinetic energy for the $n=1-3$ components during a sawtooth event in the presence of REs, considering (a) $n_\text{p}=4$ and $8$ toroidal planes and (b) $n_\text{p}=8$ and $16$ toroidal planes. It is found that varying the number of toroidal planes from 4 to 16 demonstrates that 8 planes provide sufficient accuracy in representing the $n=1$ mode.

In figure A2, we investigate the effect of (a) the RE convection velocity ($C_\text{RE}$) and (b) the time step ($\Delta t$) on the dynamics of the $n=1$ mode. By varying the $C_\text{RE}$, it is determined that $C_\text{RE} = 50 \upsilon_\text{A0}$ is sufficient for accurately simulating the nonlinear evolution of the $n=1$ mode. In addition, adjusting the time step $\Delta t$ reveals a significant influence on the linear growth rate of high $n$ modes. However, the analysis in figures 8 to 10 indicates that the predominant high $n$ and $m$ modes only enhance the loss of energetic REs beyond the $q=1$ surface during the nonlinear phase. Therefore, employing a timestep of $\Delta t=1 \tau_\text{A0}$ does not fundamentally alter the conclusions.

  A2.  The time evolutions of kinetic energy of the $n=0-3$ components for sawtooth in the presence of REs with (a) different RE convection velocity ($C_\text{RE}$) and (b) different the time step ($\Delta t$).

DownLoad: Full-Size Img PowerPoint

Appendix B.   Mode structures

In this appendix, we investigate the poloidal harmonics contributing to the growth of the $n=1$ and $n=2$ modes after $t=1800\tau_\text{A0}$, as shown in figure A1 (a). To achieve this, we conducted a Fourier spectral decomposition of the perturbed stream-function, $\delta U_{mn}$, at three different time slices. Here, the stream-function describes incompressible flow in the poloidal plane [43]. Figure A3 presents the radial mode structure of $\delta U_{mn}$ at the following time slices: (a1, b1) $t=750\tau_\text{A0}$ during the linear stage, (a2, b2) $t=1500\tau_\text{A0}$ just before the first crash in the nonlinear phase, and (a3, b3) $t=2375\tau_\text{A0}$ after the first crash. For reference, figures A3(a1) and (b1) display the radial mode structures of $\delta U_{mn}$ for the $n=1$ and $n=2$ at $t=750\tau_{\text{A0}}$, showing that the $(m,n)=(1,1)$ resistive kink mode is the dominant mode, which is consistent with the conclusions drawn in figure 2. In the nonlinear phase, a comparison of figures A3(a2) and (a3) demonstrates that the predominant harmonics for the $n = 1$ mode are $\left(m,n\right)\ =(2,1)$, $(m,n)=(3,1)$, and $(4,1)$, with the $(2,1)$ mode being particularly prominent, as indicated by the green curve in figure A3(a3). For the $n = 2$ mode, as shown in figures A3(b2) and (b3), the growing poloidal harmonics after $t = 1800\tau_\text{A0}$ include $m = 3$, $4$, and 5, with the dominant mode being $(4,2)$, as indicated by the cyan curve in figure A3(b3). These findings indicate that the high-$n$ and high-$m$ modes have a significant impact on the saturation level of the $(1,1)$ mode.

  A3.  The radial mode structures of $\delta U_{mn}$ at different time slices in the presence of REs: (a1, b1) $t=750\tau_\text{A0}$ during the linear stage, (a2, b2) $t=1500\tau_\text{A0}$ before the first crash during nonlinear stage, and (a3, b3) $t=2375\tau_\text{A0}$ after the first crash during the nonlinear stage.

DownLoad: Full-Size Img PowerPoint