Anisotropy of energetic particles and thermal transport on internal kink stability in HL-3

Achieving high pressure steady state plasma discharges is one of the major goals in advanced tokamak research. However, instabilities such as neoclassical tearing mode (NTM) often limit the maximum achievable plasma $β$⁠, and hence understanding the mechanisms that are associated with the onset of NTM is the key to avoid triggering it. One such mechanism that has been observed in many tokamak devices is that of large sawtooth crashes, which are capable of triggering NTMs which, in turn, can cause major disruptions.^1–3 To avoid large sawtooth crashes, one control approach that is often used in experiments is destabilization of sawtooth, i.e., reducing the period of sawtooth oscillations.^2–4 This can be achieved by either decreasing the plasma rotation,⁵ injecting fast ions with Neutral Beam Injection (NBI),^6,7 or utilizing Ion Cyclotron Resonance Heating.^8,9

The effectiveness of the above two methods can be understood by investigating the stability of the n = 1, m = 1 internal kink (IK) mode, where n and m are the toroidal and poloidal mode numbers, respectively. Typically, this mode becomes unstable when the safety factor at the magnetic axis, $q0$⁠, drops below 1. The IK instability is closely associated with the onset of sawtooth oscillations,^10–12 and destabilization of mode plays the fundamental role in triggering more frequent (shorter period) sawtooth oscillations. It holds for both control methods despite the fact that different destabilization mechanisms take place. In the first method,⁵ destabilization of the internal kink is achieved by reducing the gyroscopic effect associated with toroidal rotation.^13,14 In the second method as adopted in JET, the IK is destabilized by counter-passing particles in the counter-NBI heating,^4,7 resulting in shorter sawtooth periods. Note that not all NBI induced energetic particles (EPs) will play the destabilizing role. Trapped EPs¹⁵ and co-passing EPs¹⁶ can strongly stabilize the mode, leading to longer sawtooth period as has also been reported in JET.¹⁷ 

Sawtooth control is also an important topic for the new medium-sized tokamak HL-3, which is also equipped with NBI as an auxiliary heating source. In order to develop the control techniques and provide guidance to later experiments, comprehensive understanding of the stability of internal kink mode in HL-3 is essential. In this work, we consider two physics effects that have been reported to substantially modify stability of the internal kink mode. One is the kinetic effect of NBI induced EPs due to mode-particle resonances.¹⁸ Both co- and counter-NBI cases are considered, to compare the effect of passing particles with opposite circulating directions. Realistic equilibrium distribution functions of EPs are assumed, with slowing down distribution in the particle energy space and anisotropic in the particle pitch angle space, to correctly capture the effects of trapped and passing EPs on the mode stability. The other effect is anisotropic thermal transport,^19–21 which has also been found to be capable of significantly changing the stability of MHD modes^22–24 including the IK.²⁵ 

Beside individual effects of the aforementioned two factors, and their combination, have been reported on stability of the resistive wall mode²³ and the tearing mode.²⁶ This combined effect originates from the fact that anisotropic thermal transport can alter the mode-particle resonance condition by modifying the mode eigenvalue, which further affects the mode stability. However, in the case of internal kink mode, the combined effect of thermal transport and EPs remains unknown, motivating the present study, which will contribute to sawtooth control in HL-3 using NBI.

This work numerically investigates the effect of NBI induced EPs on the stability of IK in the presence of finite anisotropic thermal transport for a presentative HL-3 plasma. We employ the MHD-kinetic hybrid code MARS-K,²⁷ where an analytic model of the EPs equilibrium distribution was implemented.¹⁸ A non-perturbative approach is adopted in MARS-K to self-consistently determine the modification to the mode eigenvalue as well as the mode structure. Both adiabatic and non-adiabatic contributions of EPs will be considered, for both passing and trapped EPs though the latter has small fraction with the tangential NBI. Here, the non-adiabatic effect of EPs is due to resonances between the mode and different drift frequencies of EPs (essentially wave–particle Landau resonances). In this work, the transit, bounce, and precessional drift resonances will all be considered.

The rest of the paper is structured as follows. Section II briefly introduces the MARS-K computational model, including implementation of the EP distribution function and the anisotropic thermal transport model. The HL-3 equilibrium considered in this work is also introduced in this section. Results from toroidal computations with MARS-K are reported in Sec. III. The conclusion is drawn in Sec. IV.

Various forms of distribution can be specified by adopting different sets of coefficients $Ci, ζi$ in Eqs. (8) and (9). Table I lists the $Ci, ζi$ assumed in this work to describe EPs from co- and countercurrent tangential NBIs. The EP distribution function is shown in Fig. 1(a), assuming $δζ02=0.1$⁠, $ζ0$= 0.64 (co-NBI) and −0.64 (counter-NBI), to match the predicted distribution from the TRANSP simulation performed for HL-3.²⁹ Shown in Fig. 1 are also the radial profiles of the surface-averaged EP pressure and number density, defined by the ratios $PEP/Pth$ and $NEP/Nth$⁠, respectively. Here, the subscript “th” denotes thermal particles.

Note that the perturbed distribution function used to calculate $p∥EP$ and $p⊥EP$ also depends on the perturbed fluid quantities such as the plasma displacement and the perturbed magnetic field. This consequently forms a closed MHD-kinetic hybrid model which self-consistently determines the mode eigenvalue and structure.

We have discussed anisotropic aspects associated with the EPs. One is the equilibrium distribution of EPs in the particle phase space, which is anisotropic due to the nature of the NBI. The other is the perturbed pressure tensor due to EPs (⁠ $p∥EP$ and $p⊥EP$⁠), which is anisotropic as a result of the drift kinetic calculation. On top of these, the present study also considers another type of anisotropy, associated with thermal particles. The latter appears in the last two terms from the right-hand side of Eq. (5), i.e., the anisotropic thermal transport terms (TTTs). The anisotropy of thermal transport is considered here to be more relevant to realistic tokamak conditions, where the parallel thermal conduction coefficient, $χ∥,$ is several orders of magnitude larger than that in the perpendicular direction, $χ⊥$⁠. For most of modeling results presented in this work, we fix the ratio $χ∥/χ⊥=106$⁠, with $χ⊥=3×10−6$ and $χ∥=3×100$ normalized by $a2/τA$⁠. Here, $a$ is the minor radius of the plasma and $τA$ the toroidal Alfvén time. Later in this study, this ratio will also be varied to investigate its impact on the stability of the internal kink mode.

The plasma equilibrium considered in this work is obtained by the TRANSP simulation, with the plasma cross section shown in Fig. 2(a). This HL-3 plasma has major radius $R0=1.8 m$⁠, aspect ratio $R0/a=3$⁠, and the toroidal magnetic field $B0=1 T$⁠. The equilibrium pressure and safety factor profiles (with the on-axis safety factor $q0=0.8)$ are plotted in Figs. 2(b) and 2(c), respectively. The on-axis safety factor will be scanned later on with the corresponding radial profile self-consistently determined using a fixed boundary equilibrium solver (the CHEASE code). Subsonic plasma toroidal flow is also considered in this study, with the radial profile of the angular frequency shown in Fig. 2(d).

Before showing the modeling results, we make a remark on the relation between the parallel thermal transport and the plasma equilibrium. Normally, certain equilibrium quantities (e.g., the plasma pressure and safety factor as shown in Fig. 2, as well as the plasma density and temperature) are functions of magnetic flux, which follows from long-time equilibration as particles move along the magnetic field lines. For thermal particles, this parallel motion is characterized by the sound speed and thus the parallel thermal transport is not infinite, but it can indeed be very large as we assume in this work. Other quantities, such as the plasma current density or poloidal magnetic field, are not functions of the flux surface even in the equilibrium state. All perturbed quantities in our model are not functions of the equilibrium flux surface. The anisotropic transport terms are only associated with these perturbed quantities in our linearized model.

We start by examining the effect of EPs on the stability of the internal kink mode including anisotropic thermal transport. The MARS-K results are reported in Fig. 3, where variations of the internal kink mode growth rate with increasing on-axis safety factor are compared among the fluid, co- and counter-NBI cases, as well as between cases with and without thermal transport. Note that the IK instability is primarily caused by the fact that the on-axis safety factor $q0<1$ (i.e., current-driven), as shown by the black line (labeled as “fluid”) in Fig. 3(a) where neither the kinetic effect of EPs nor the thermal transport is present. On the other hand, kinetic effects of EPs and the thermal transport do significantly affect the mode stability. These additional effects beyond the standard fluid model are the focus of the present study.

In the absence of thermal transport, compared to the fluid model, the effects of EPs in both co- and counter-tangential NBI are similar—they both play a stabilizing role. However, with inclusion of the thermal transport effect, the roles become opposite as demonstrated in Fig. 3(b). EPs induced by co-NBI still stabilize the IK for $q0$ < 0.95. On the other hand, EPs with counter-NBI case strongly destabilize IK. Note that this combined effect of thermal transport and EPs on the mode is self-consistently determined by simultaneously including the two effects in the MARS-K model.

Both adiabatic and non-adiabatic effects of EPs are included as follows:

1. adiabatic contribution from both passing and trapped EPs, which will be referred as “AP” and “AT,” respectively;

2. non-adiabatic contribution from passing EPs, which is due to the transit resonance with the mode (referred as “NP”);

3. non-adiabatic contribution from trapped EPs, due to both bounce (“NTB”) and precession (“NTD”) drift resonance.

All these portions contribute to the total kinetic effect of EPs reported in Fig. 3. We wish to understand which one of the above is responsible for the opposite effect in the two NBI cases as shown in Fig. 3(b). This is feasible since MARS-K allows including the individual contribution from each part. The result of such detailed modeling is shown in Fig. 4, where variations of the growth rate and real frequency of the mode are plotted, when the aforementioned kinetic contributions are progressively included for both the co- and counter-NBI cases, following the order “AP $→$ AT $→$ NP $→$ NTB $→$ NTD,” as explicitly shown in Table II.

For the case without thermal transport, the adiabatic contribution of passing EPs (AP) is stabilizing, while that of trapped EPs (AT) is destabilizing. Non-adiabatic contribution from passing (NP) and trapped EPs (NTB and NTD) has weak effect on the mode stability. More importantly, different kinetic contributions have similar effects in both co- and counter-NBI cases, leading to the similar stabilizing effect in Fig. 3(a). When thermal transport in included, although the kinetic contributions from other parts remain similar between the two NBI cases, the non-adiabatic contribution from passing EPs due to transit resonances becomes opposite between the two NBI cases, i.e., being strongly stabilizing with the co-NBI cases but destabilizing with the counter-NBI case. This explains the different effects on mode stability between the two NBI cases as shown in Fig. 3(b). This is reasonable, since the large difference cannot come from trapped EPs because they have the same distribution in the two NBI cases as demonstrated in Fig. 1(a), resulting in the same stabilizing contribution.

Equation (10) appears to show that inclusion of the thermal transport effect does not directly affect the resonance condition. However, it can modify the condition in an indirect manner, i.e., by altering the mode eigenvalue. On the other hand, the circulating direction of passing particles also affects the resonance condition, which will further contribute to the different non-adiabatic contribution of passing EPs with the two NBI cases.

Figure 4(b) shows the real frequency of the mode, for the same four cases as in Fig. 4(a). In general, the non-adiabatic contribution of trapped EPs tends to increase the mode frequency, while that of passing EPs decreases it. The latter is more pronounced in the counter-NBI case in the presence of thermal transport, leading to lower frequency in this case when all kinetic effects are included.

The difference between the IK stability with the co- and counter-NBI cases also depends on the ratio $χ∥/χ⊥$⁠, i.e., the anisotropy of the thermal transport. This is demonstrated in Fig. 5, where variations of mode eigenvalue are plotted while scanning the mentioned ratio at two different $χ⊥$⁠. Despite the fact that the growth rate is higher with larger $χ⊥$⁠, we notice similar variations of mode eigenvalue with respect to $χ∥/χ⊥$ at both $χ⊥$ values. The growth rate is reduced when decreasing ratio $χ∥/χ⊥$ in both NBI cases and becomes similar at small $χ∥/χ⊥$⁠. This trend implies that in the limit of vanishing anisotropy, i.e., isotropic thermal transport where $χ∥/χ⊥=1$⁠, the kinetic effects of EPs will provide similar stabilizing effect in both the co- and counter-NBI cases. With increasing ratio, due to its indirect modification to the resonant condition via mode eigenvalue, the mode becomes more unstable, in particular with the counter-NBI case. Consequently, the difference between the kinetic effects of EPs with the two NBI cases becomes stronger at larger $χ∥/χ⊥$⁠.

Since we adopt the non-perturbative method in MARS-K, we are also able to study the modification to mode structure due to inclusion of the thermal transport effect. The results are reported in Fig. 6, where radial profiles of the plasma radial displacement are compared with and without thermal transport and with the co- and counter-NBI cases, respectively. Modification of the mode structure with inclusion of thermal transport is weak for the co-NBI case, but it is strong for the counter-NBI case, especially inside the q = 1 rational surface where the mode amplitude is relatively smaller in the presence of thermal transport. Note that this modification is mainly contributed by the non-adiabatic contribution from passing EPs.

We have investigated effects of the NBI induced EPs on the linear stability of internal kink mode in the presence of anisotropic thermal transport. In the absence of the thermal transport effect, MARS-K modeling shows that the kinetic effect of EPs with the co- and counter-NBI cases has similar effects on the internal mode, i.e., both being stabilizing. However, inclusion of anisotropic thermal transport results in opposite kinetic effects of EPs on the mode stability between the two NBI cases: stabilizing with the co-NBI case and destabilizing with the counter-NBI case. These findings qualitatively agree with observations in the JET experiments.¹² 

Detailed investigation reveals that the large difference mainly comes from the opposite non-adiabatic contribution of passing EPs due to transit resonances. This can be understood by examining the resonant operator, where both the circulating direction of passing EPs and mode eigenvalue matter. The anisotropic thermal transport can indirectly modify the resonance condition by significantly altering the mode eigenvalue, which further contributes to the different kinetic effect of EPs with the two NBI cases. It is therefore expected that varying the ratio between the parallel and perpendicular thermal transport coefficients can alter the kinetic effect of EPs. This is indeed the case, as the MARS-K computations show that the differences between kinetic effects of EPs between the two NBI cases are negligible for smaller thermal anisotropy ratio, but become significant when increasing the ratio (i.e., the anisotropy of thermal transport). For realistic ratio values which (> $106$⁠) in tokamak plasmas, different kinetic effects of EPs between the two NBI cases can be important. Similar differences between the two NBI cases, in the presence of thermal transport, are also observed in the computed mode structure, where modifications to the mode structure with inclusion of thermal transport are more pronounced with the counter-NBI case.

The conclusion drawn here can have further implications for sawtooth control in later HL-3 experiments. It is expected that the sawtooth period would be lengthened for the on-axis co-NBI due to the stabilizing effects of EPs on the internal kink mode and be shortened in the counter-NBI case. This in turn indicates that the counter-NBI may be more beneficial in terms of producing frequent sawtooth (without considering the prompt EP loss issue).

As a final remark, the results reported here are limited in the following aspects: 1) Due to linear nature of MARS-K, we only consider the linear stability of internal kink mode, without nonlinear growth and saturation of the mode. 2) Only the linear effect of thermal transport, i.e., competing with compression terms in Eq. (5) to determine perturbed pressure, is considered, while its nonlinear effects such as local flattening of the pressure profile are not captured. 3) We have ignored the finite orbit correction effect for EPs, which can be particularly important for passing particles as long as the internal kink instability is concerned.³⁰ 4) We have focused on the kinetic effect of EPs on the stability of internal kink mode. On the other hand, internal kink mode induced 3D perturbations may also impact the distribution and confinement of EPs, again making the problem nonlinear, which is beyond the scope of the present investigation.

This work was supported by the National Natural Science Foundation of China under Grant No. 12105082.

The authors have no conflicts to disclose.

X. Bai: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Resources (lead); Visualization (lead); Writing – original draft (lead). Y. Liu: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). G. Dong: Data curation (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal). G. Z. Hao: Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal).

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