Magnetic island formation and rotation braking induced by low-Z impurity penetration in an EAST plasma

Impurity injection introduces two direct modifications to the equilibrium profiles (figure 3). On the one hand, the impurity radiation cooling leads to the temperature profile contraction; on the other hand, the impurity ionization increases the electron density at the impurity cold front, which is identified from the local peak of the electron density profile (figure 3(a)). As a consequence, the pressure profile steepens and its gradient is enhanced at the impurity cold front around the q = 2 surface (figure 3(b)). This local variation in pressure gradient is predominantly from the temperature profile change in the simulation case considered here.

Zoom In Zoom Out Reset image size

Despite the dynamic nature of the impurity injection process, the profile evolution is quasi-static and the force balance ${\mathrm{d}}p/{\mathrm{d}}r \approx J \times B$ is well maintained (figure 4). The corresponding n = 0 component of the parallel current perturbation $\delta J_{\parallel}$ mainly comes from the perturbed Pfirsch–Schlüter (PS) current $\delta J_{\mathrm{ps}} = -2 \frac{1}{B_{\mathrm{p}}} \frac{r}{R} \frac{{\mathrm{d}}\delta p}{{\mathrm{d}}r} \cos\theta$ in the outer region on both sides of the q = 2 surface, which is induced by the pressure gradient perturbation ${\mathrm{d}}\delta p/{\mathrm{d}}r$. However, the parallel current perturbation peaks around the q = 2 surface can not be accounted by the perturbed Pfirsch–Schlüter current $\delta J_{\mathrm{ps}}$ alone.

Zoom In Zoom Out Reset image size

The local positive current perturbation peak denoted by the red arrow correlates well with the local impurity radiation peak (figure 5). Similar skin current structure due to edge radiative cooling is also reported in previous studies [21–23]. The local current perturbation dip denoted by the blue arrow collocates with the surge of plasma resistivity in the cold region and is dominated by the current diffusion (figure 5). It is important to note that the location of impurity radiation peak is usually different from that of the enhanced plasma resistivity, because the impurity radiation is dominated by the line radiation, which is the strongest around $T_{\mathrm{e}} \sim 10^2\,\mathrm{eV}$, whereas the resistivity is proportional to $T_{\mathrm{e}}^{-3/2}$ and more enhanced in the colder region.

Zoom In Zoom Out Reset image size

The TM growth are driven by the parallel current perturbation as expected and confirmed from the correlation observed for the $2/1$ as well as the $3/1$ modes (figures 7(a) and (b)). Based on previous developed analytic theories on the nonlinear resistive growth of magnetic island including the interchange effects, the following modified Rutherford equation (MRE) may be introduced for discussion [24–27]

$$\begin{equation} \frac{1}{\eta} \frac{\partial w}{\partial t} \sim \Delta^{^{\prime}} + \frac{\alpha_1 D_{\mathrm{R}} + \alpha_2 E_{\phi_0}\delta \eta_{m,n}}{\sqrt{w^2 +w_d^2}} \end{equation} \tag{ 7 }$$

where $\Delta^{^{\prime}}$ is the tearing instability index, $\delta \eta_{m,n}$ is the (m, n) helical component of the perturbed resistivity, and $E_{\phi_0}$ is the toroidal electric field due to loop voltage. Here $D_{\mathrm{R}} \approx \frac{\epsilon^2_s\beta_p}{s} \frac{L_q}{L_p} \left( 1-\frac{1}{q^2} \right)$ is the resistive interchange parameter with $\epsilon_s = r_s/R_0, \beta_p = 2\mu_0p/B_p^2, L_q = q/q^{^{\prime}} = r_s/s$, and $L_p = p/p^{^{\prime}}$, which is usually negative in a tokamak with monotonic safety factor q profile [24], $w_d = 2\sqrt{2}\left( \kappa_{\perp}/\kappa_{\parallel} \right)^{1/4} \left(r_sR/ns \right)^{1/2}$ is the finite thermal diffusion length scale [28], and the coefficients α₁ and α₂ may be determined from the more quantitative theory or fit from the corresponding simulation or experimental results. The overall time evolution of all key terms in equation (7) are summarized in figure 6, which we delineate further below in details.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the $2/1$ mode, the surface-averaged parallel current perturbation around the q = 2 rational surface $\left\langle \delta J_{\parallel} \right\rangle_{q = 2}$ increases rapidly at $t = 1.5\,\mathrm{ms}$, and the $2/1$ island width grows up thereafter at $t = 2.5\,\mathrm{ms}$ (figure 7(a)). The parallel current perturbation $\left\langle\delta J_{\parallel} \right\rangle_{q = 2}$ is primarily induced by the enhanced impurity radiation, as indicated by its correlation with the local radiation power in figure 8(a). The $m = 0,n = 0$ component of the parallel current density perturbation $\left\langle \delta J_{\parallel}\right\rangle$ drives the tearing growth through the $\Delta^{^{\prime}}$ in the Rutherford equation, and the linear TM growth is just the results of major change to the current profile, i.e. $\left\langle \delta J_{\parallel}\right\rangle$, which itself is induced by the impurity radiation. The helical current terms in the MRE (7) contribute to the nonlinear island growth, but not the trigger to the initial TM onset. In absence of any impurity deposition, the equilibrium adopted in the simulation is stable to the linear TM, i.e. $\Delta^{^{\prime}}\lt0$. Only after the impurity deposition is introduced to the equilibrium and the subsequent penetration arrives at the corresponding rational surfaces, the spontaneous TM island growth are mainly driven by the now unstable parallel current profile with $\Delta^{^{\prime}}\gt0$ as direct results of impurity radiative cooling. It is the impurity radiation effects on the current profile that triggers the linear growth of a TM. Once into the nonlinear stage, the thermal-resistive mechanism involving power balance between Ohmic heating and impurity radiations starts to contribute to the island development. The nonlinear island growth in our simulations here is not of neoclassical tearing mode (NTM) but the resistive tearing, since neither the neoclassical closure nor the bootstrap current is included in these simulations.

Zoom In Zoom Out Reset image size

After $t = 4.5\,\mathrm{ms}$, the current perturbation becomes dominated by resistive diffusion due to the enhanced plasma resistivity. The $1.0 \,\mathrm{ms}$ delay (t = 1.5–2.5 ms) in the mode growth may be due to the well-known stabilization effect from flux-averaged pressure gradient and magnetic curvature [24, 29], which reduces the cylindrical tearing instability parameter $\Delta^{^{\prime}}$ to an effective value [26], i.e. $\Delta_{\mathrm{eff}}^{^{\prime}} = \Delta^{^{\prime}} + \sqrt{2}\pi^{3/2} D_{\mathrm{R}}/w_d$. This derives from the fact that perturbed pressure leads to a parallel current perturbation outside the resistive layer through magnetic curvature, which contributes to the jump in the logarithmic derivative, i.e. the $\Delta^{^{\prime}}$. From figure 8(a), it is clear that the local pressure gradient is enhanced around $1.5\,\mathrm{ms}$ and decreases toward zero right before the mode growth ($t = 2.5\,\mathrm{ms}$). More importantly, the absolute value $|D_{\mathrm{R}}|$ indeed increases to a larger value right before the mode growth, which represents the stabilization effect, and drops rapidly once the mode begins to grow (figure 8(b)). With the impurity inward penetration, the current perturbations move along with the cold front and cross different rational surfaces to trigger TMs with multiple helicities.

A stationary equilibrium plasma framework is adopted in the simulations, i.e. $V_0 = 0\,\mathrm{m}\,\mathrm{s}^{-1}$, to study the impurity effect on the rotation of magnetic islands alone, and the plasma equilibrium dose not start rotating without the impurity injection. A small perturbed toroidal velocity at the beginning of simulation is caused by the impurity injection, then it increases rapidly to its peak and decays slowly thereafter (figure 9). The mode frequency of the n = 1 normal component of the perturbed magnetic field $B_{\mathrm{r}}$ is slowly damped from beginning, and the magnetic islands become almost stationary after $t \gtrsim 4\,\mathrm{ms}$ while the toroidal flow still remaining finite. The toroidal rotation frequency on the q = 2 surface is approximate twice that on the q = 3 surface, whereas the mode frequencies on these two surfaces are almost same during t = 0–4.0 ms before locking. Such a mode excitation with finite frequency and the subsequent gradual damping and eventual locking is also observed in experiments, however, such a mode locking is previously attributed to the electromagnetic torque braking due to the error field from the tungsten protector limiter [14]. The local magnetic perturbation on the q = 2 surface shows the initial dominant poloidal mode number m = 3 before $t = 2.5\,\mathrm{ms}$, which is later replaced by the m = 2 poloidal component only after the excitation of the $2/1$ mode by the arrival of the impurity radiation peak on the q = 2 surface (figure 9(a)).

Zoom In Zoom Out Reset image size

In EAST experiments, the decrease in rotation velocity is attributed to the mechanism of mode locking following the redistribution of low-Z impurity concentration, which primarily contributes to the increase in the electromagnetic brake torque, through the increased island size by impurity radiation cooling as well as the reduction of rotation velocity itself [14]. Our simulation data indicates that the rotation in the subsequent simulation is induced first by the poloidal and toroidal pressure gradient as a direct result from the localized impurity radiative cooling. Once the impurity penetration reaches the rational surface where a TM island is driven to grow, the electromagnetic torque becomes the dominant driver for the island region plasma rotation, which can be seen from the time history of the poloidal components of the velocity, the pressure gradient, and the electromagnetic torque averaged over the q = 2 surface (figure 10). The electromagnetic torque includes contributions from both the diamagnetic current due to the change in radial pressure gradient and the inductive current due to the TM growth. The enhanced diamagnetic drift is equivalent to the enhanced diamagnetic current component of the electromagnetic torque, which contributes to the rotation change of the island region plasma rotation. The rotation eventually slows down after the impurity penetration induced pressure gradient and electromagnetic torque gradually decrease and vanish.

Zoom In Zoom Out Reset image size

The impurity radial transport mainly depends on diffusion, and its penetration in the poloidal direction is dominated by the perturbed poloidal flow through convection. Meanwhile, the localization of the impurity density along with the radiation power not only lead to the island formation, but also correlates to the phase of the TM. As indicated by figure 11(c), the localization of the impurity radiation power corresponds to the O-point of the $2/1$ island, which is widely observed both in experiments and simulations. The radiation leads to temperature cooling that enhances the plasma resistivity, thus reduces the current density and causes the O-point of the island. And the island may confine the impurity ions locally due to its magnetic structure. Therefore, the magnetic island rotates along with the impurity density in the poloidal plane, which is influenced by the perturbed poloidal flow, as identified by the correlation between the poloidal distribution of the radiation power, the perturbed poloidal flow, and the n = 1 mode phase evolution (figures 11(a) and (b)).

Zoom In Zoom Out Reset image size

The island rotation amplitude can be measured by the integral of the perpendicular vortex associated with the mode in the poloidal plane, which decreases toward zero gradually (figure 12(a)). Whereas the island rotation velocity is a three-component vector with a 3D distribution in space, the integral of the perpendicular vortex $\int U_{\perp} {\mathrm{d}}S = \int \nabla_{\perp} V {\mathrm{d}}S$ is a single-valued global measure of the total degree of island rotation. The definition or the measurement of such a quantity does not rely on the presence or absence of a mode in the velocity field. Meanwhile, the impurity penetration front can be indicated by the local enhanced electron density peak, which sweeps inward in a step-wise manner across rational surfaces. The mode rotation amplitude rapidly shoots to its peak value in the beginning when the impurity is localized in the bottom region with a strong up-down asymmetric distribution, then gradually slows down as the impurity penetrates inward along with more uniform toroidal and poloidal spreading. This simulation results agree the experimental observation that the modes can be locked following the redistribution of the low-Z impurity concentration [14], despite the fact that there is no error or external magnetic field at the perfectly conducting wall boundary in simulations. Similarly, Joint TEXT tokamak (J-TEXT) experiments demonstrate strong correlation between the $2/1$ TM rotation and impurity distribution as well [30, 31].

Zoom In Zoom Out Reset image size

The inclusion of higher-n modes in simulation accelerates the mode rotation drop toward zero and the subsequent stationary state in the simulation (figure 12(b)). This suggests that the high-n helical structures may be able to introduce additional braking effects, likely through the electromagnetic torques from the magnetic island chains, as well as the coupling and overlapping of magnetic islands on the neighbouring rational surfaces.

Higher-n modes also significantly impede the impurity inward penetration and as a consequence the radiation peak stays longer upon the q = 2 surface (figure 12(c)), which could be due to a combined effects from the stochastic field [32] and the $2/1$ magnetic island itself [33]. By contrast, in the simulation case with the n = 0 component only, the impurity front almost directly penetrates into the central region in absence of magnetic islands. This agrees with the observations that during an MGI experiment, the impurity penetration usually stops along the q = 2 surface [20, 34], and such an agreement highlights the critical roles of the higher-n modes in the impurity penetration process.