A new mechanism for internal kink stabilization by plasma flow and anisotropic thermal transport

Stability of the internal kink (IK) mode is numerically investigated with inclusion of the anisotropic thermal transport effect and the toroidal plasma flow. It is found that anisotropic thermal transport, in combination with plasma flow, stabilizes the IK in two ways. One is direct stabilization of the mode synergistically with plasma flow. The other, indirect stabilization involves generation of a finite mode frequency in static plasmas by thermal transport, which in turn invokes wave-wave resonance damping of the mode via interaction with stable shear Alfvén waves. This second IK stabilization mechanism is further corroborated by examining the eigenmode structure, which peaks at the radial locations where the mode frequency matches that of the shear Alfvén wave. Finally, two branches of unstable IK are identified, with mode conversion occurring at certain plasma flow speed and thermal transport level. These findings provide new physics insights in the IK stability in tokamak fusion plasmas.

(Some figures may appear in colour only in the online journal) Sawtooth oscillation is one of the universal phenomena that has been observed in many tokamak devices, where periodic collapse of the plasma core density and temperature occurs [1].Though frequent and small sawtooth can be beneficial for preventing impurity accumulation in the plasma core [2], large ones such as the so-called monster sawtooth can be dangerous to the plasma operation by triggering the Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.neoclassical tearing mode (NTM), which degrades the plasma confinement and may eventually lead to plasma disruption [3,4].To avoid seeding NTMs, it is therefore of critical importance to understand the sawtooth oscillation physics and the magnetohydrodynamic (MHD) instability underlying these oscillations [5][6][7], as well as to implement control techniques to avoid large-amplitude sawtooth crashes [8].Extensive theoretical and experimental efforts have been devoted for these purposes.It is now widely accepted that sawtooth is triggered by the unstable m/n = 1/1 internal kink (IK) component [9] (where m and n are, respectively, the poloidal and toroidal mode numbers), whose stability has been proven to be crucial in explaining the sawtooth behavior in experiments.
The n = 1 IK mode is the most important MHD instability in the plasma core region.Typically, the mode becomes unstable when the safety factor at the magnetic axis, q 0 , drops below unity.Stability of the IK has been extensively investigated in different magnetic geometries and under various plasma conditions [9][10][11].Many physics effects have been found to significantly affect the IK stability, including kinetic effects of energetic particles (EPs) [6,12,13] and gyroscopic stabilization due to toroidal rotation of the plasma [14,15].Energetic ions, such as fusion born alpha particles and trapped fast ions produced by auxiliary heating, can stabilize IK, while counter passing EPs can destabilize it as observed in JET experiments [16].Toroidal plasma rotation generally provides a stabilizing effect on the mode.The flow shear is also found to play an important role when combined with drift kinetic stabilization of trapped EPs [15].
The conventional understanding of the IK stabilization by the plasma toroidal rotation relies on the gyroscopic effect, which requires a rather fast flow (close to the sound speed) for the effect to become prominent.This rotational stabilization mechanism is very different from that for other MHD instabilities such as the resistive wall mode (RWM) [17], where strong continuum wave resonance damping occurs at relatively slow plasma flow (well below sound speed).The ultimate reason for the strong continuum wave damping (on top of the potential drift kinetic damping from plasma particle species) of the RWM is the fact that the mode is essentially locked to the resistive wall and thus does not rotate together with the plasma, creating resonance conditions with (stable) continuum waves in the plasma.The IK, being a core localized instability that weakly couples to the outer wall, on the other hand generally rotates together with the plasma (i.e.being static in the plasma frame).This leads to difficulty to achieve resonance damping of the IK by coupling to stable continuum waves of the plasma.
However, if certain additional physics effect can alter the IK mode frequency in the plasma frame, resonant damping with plasma continuum waves is again possible, which in turn modifies the mode stability.One example of such an additional physics is the kinetic effect of EPs, which can trigger the fishbone instability with finite frequency even in a static plasma.As a consequence, continuum wave resonance can occur for this kind of finite-frequency mode [18].This work reports a new physics mechanism of generating finite mode frequency for the IK (even in a static plasma), which is in turn subject to shear Alfvén wave continuum damping.This physics is associated with the anisotropic thermal transport, which was also previously found to alter the mode frequency for the RWM [19].
We consider linear stability of the resistive IK in toroidal geometry, taking into account both toroidal flow and the anisotropic thermal transport effects.We utilize the MARS-F code [20] to numerically compute the mode stability.We find that thermal transport can stabilize the IK in two ways: one is to directly stabilize the mode together with toroidal flow; the other is an indirect way via producing a finite mode frequency in the plasma frame, and thus allowing the mode stabilization due to resonances with shear Alfvén waves.The second mechanism is the key result of the present work.This new stabilizing mechanism, together with other stabilizing physics such as that identified by Bussac [10], may offer a better explanation of the experimental observation that the on-axis safety factor can be robustly below unity in certain DIII-D plasmas [21].
We consider the effect of anisotropic thermal transport on the stability of IK by utilizing the MARS-F code that solves the full toroidal, single fluid, linearized MHD equations including the equilibrium toroidal flow and the anisotropic thermal transport terms (TTT) for the perturbed pressure equation ρ + ∇ • (ρ eq ξ) = 0 (1) where γ is the (generally complex) eigenvalue of the mode.inΩ represents the Doppler shift due to toroidal plasma rotation.R is the major radius of the torus, Ẑ and φ are the unit vectors along the vertical and toroidal directions, respectively.The plasma displacement, perturbed velocity, magnetic field, density and plasma current are denoted by ξ , v, b, ρ, j, respectively, whilst ρ eq , B eq , J eq and P eq represent the equilibrium plasma density, magnetic field, current and thermal pressure.
The anisotropic thermal transport enters the above model via the perturbed pressure equation ( 5) as the last two terms from the right-hand side of the equation.Note that these terms are absent in standard single-fluid MHD equations.These terms have been implemented in MARS-F following the Braginskii closure [22] but with the transport coefficients directly applied to the perturbed pressure in the single fluid model, as has also been adopted in previous work [23,24].
We adopt the anisotropic transport model due to the fact that the parallel and perpendicular thermal transport coefficients χ ∥ and χ ⊥ are drastically different in a typical tokamak plasma.The ratio between the parallel and perpendicular thermal conduction coefficient, χ ∥ /χ ⊥ , is expected to be 10 6 ∼ 10 10 in a typical tokamak.In this work, we will assume χ ∥ /χ ⊥ = 10 6 in most cases.Apart from their ratio, the values of χ ∥ and χ ⊥ also matter.Reported in this work are the normalized values for these coefficients, by a 2 /τ A with a being the plasma minor radius and τ A the toroidal Alfvén time.In this study, these thermal transport coefficients are assumed constants that do not vary along the plasma radial coordinate.
We note that the most well-known effect of anisotropic thermal transport in tokamak plasmas is to produce local flattening of the pressure profile inside the magnetic island.This is a non-linear effect that is not considered in the present study.The linearized TTT also affect the linear stability of MHD modes by competing with the sound wave propagation [24]this is the effect we capture in this work.
The plasma considered in this work is from the HL-3 tokamak, with cross section shown in figure 1(a).The plasma major radius, aspect ratio and the toroidal magnetic field are R 0 = 1.80 m, R 0 /a = 3, B 0 = 1 T, respectively.The equilibrium pressure and safety factor profiles are plotted in figures 1(b) and (c) for two representative cases, with the onaxis safety factor values of q 0 = 0.9 and 0.97.Later on, we will scan q 0 and obtain the self-consistent toroidal equilibrium for each q 0 value.The radial profile of the toroidal rotation of the plasma is shown in figure 1(d), with the overall amplitude again being scanned.Note that this study assumes static (Grad-Shafranov) equilibrium.The effect of the toroidal plasma flow is only reflected in the perturbed MHD equations.It is well-known that toroidal flow also modifies the plasma equilibrium by a factor of ∼exp(M 2 ) with M being the sound Mach number [25].In our work, we consider toroidal flow speed that is well sub-sonic (M < 0.05), meaning that the flow correction to the equilibrium is about 0.25% which we deem as negligible.
We start by reporting the first stabilization mechanism for the n = 1 IK by anisotropic thermal transport.Figure 2 compares the MARS-F computed mode growth rate while varying the on-axis safety factor, between two cases-one without and one with finite plasma toroidal flow.The thermal transport effects are similar in both cases-the IK is destabilized when q 0 is well below unity but the mode stabilization starts to take place as q 0 exceeds a critical value.This stabilizing effect is more significant with higher level of thermal transport.The likely physics mechanism behind this stabilization is the following.The stability of the IK is sensitive to the pressure gradient and any physics effects that alter the pressure gradient near the q = 1 rational surface will have significant effect on the mode stability.It is also confirmed by MARS-F modelling that the growth rate of the mode is reduced if the pressure profile is locally flattened near the q = 1 surface.In linear theory, even though the anisotropic thermal transport does not directly flatten the pressure profile (as what would happen in nonlinear theory), it does modify the stability drive due to the pressure gradient by competing the sound wave propagation which in turn determines the perturbed pressure.This produces a stabilizing effect on the IK.
Despite the similarity of the results reported in figures 2(a) and (b), one striking difference is the complete suppression of the IK instability at q 0 > 0.97 as shown in figure 2(b), with finite plasma flow and sufficiently high thermal transport.This complete stabilization is achieved by the stabilizing effect of anisotropic thermal transport synergistically with plasma flow.) and the on-axis safety factor ( q 0 = 0.97).Shown are two branches of the mode, indicated by red-solid and blue-dashed lines, respectively.
To better understand this synergy effect, we also scan the toroidal rotation frequency at fixed q 0 = 0.97 (figure 3).Two interesting observations are made here.First, MARS-F finds two unstable IK branches with varying rotation frequency.(MARS-F employs an eigenvalue solver that allows finding multiple branches of the instability upon proper initial guess.).The first unstable branch (branch 1 shown in red) is eventually stabilized with increasing rotation.This is the same solution reported in figure 2, with the rotation threshold for achieving full stabilization at ω 0 /ω A = 0.0085 when χ ⊥ = 10 −4 .Note this threshold value is substantially subsonic.The stabilization is therefore not associated with the gyroscopic effect but rather with the thermal transport effect.The real frequency of this branch first increases with the plasma rotation frequency and then decreases, reaching a peak value around ω 0 /ω A = 0.005.
The second, even more interesting, observation from figure 3 is the other unstable branch (branch 2 shown in blue).This branch is destabilized with increasing rotation.More importantly, complete stabilization of this branch occurs at slow plasma flow (ω 0 /ω A < 0.0041 in our case).Evidently, this stabilization cannot be explained by the synergetic effect of thermal transport and plasma flow, for it only occurs at vanishing or slow plasma flow.On the other hand, we find that with branch 2, anisotropic thermal transport induces finite mode frequency even in a static plasma.This finite mode frequency (in the plasma frame) then creates resonance conditions with plasma continuum waves, such as the sound and shear Alfvén waves.Such wave-wave interactions with stable continua can then facilitate the IK mode stabilization.Although the mode frequency are significantly different for the two branches, they both peak around ω 0 /ω A = 0.005.As will be shown later, this is the region where the mode conversion occurs.
This type of resonance is in fact confirmed by carefully examining the m/n = 1/1 eigenmode structure (figure 4).Two representative cases are shown in figure 4: (a 3) and (b) vanishing plasma flow (corresponding to branch 2 from figure 3).The latter is also plotted along the safety factor q in the lower panel of (c), together with the shear Alfvén wave frequency k ∥ v A (upper panel) near the q = 1 surface (indicated by vertical dashed lines in all plots).The horizontal line in the upper panel of (c) indicates the computed mode frequency in the plasma frame.The two solid vertical lines in (c) indicate the radial locations where the mode frequency matches the shear Alfvén wave frequency.Assumed is a case with finite thermal transport coefficients (χ ⊥ = 10 −4 , χ ∥ /χ ⊥ = 10 6 ) and an equilibrium with the on-axis safety factor of q 0 = 0.97.evident that different mechanisms are involved in the mode stability for branches 1 and 2. While branch 1 represents a typical IK structure (figure 4(a)), branch 2 is substantially different (figure 4(b)).In particular, the plasma radial displacement for branch 2 has double-peaking around the q = 1 surface, resulted from resonances with plasma waves.Note that since no EPs are included in the present modelling, drift resonances with EPs [18] are therefore excluded.The nature of the resonance is further revealed by comparing the radial location of the double-peak of the plasma displacement with the frequency matching location between the mode and the shear Alfvén wave (figure 4(c)).The upper panel of figure 4(c) plots the radial distribution of the shear Alfvén wave frequency ω A near the q = 1 surface.The resonance condition follows where v A is the Alfvén speed evaluated at the magnetic axis and k ∥ the parallel wave number.The resonance occurs when the mode frequency ω r (defined in the plasma frame), indicated by the horizontal line in the upper panel of figure 4(c), matches ω A .Note that we have also checked the sound wave frequency but the latter is too slow to match the mode frequency in our case.We emphasize that the finite mode frequency ω r here is induced solely by the thermal transport effect.Both the upper and lower panels of figure 4(c) show clear alignment of the shear Alfvén wave resonance locations and the double-peak locations in the computed IK eigenmode.This is the direct evidence that these eigenmode peaks are the consequence of the mode resonance with the shear Alfvén wave.It is also this resonance that leads to complete stabilization of the IK (i.e.branch 2) at vanishing plasma flow.We also plot the eigenvalue of both IK branches in the complex plane (figure 5) while varying both the plasma rotation frequency and the thermal transport coefficients.It is evident that these two branches are located in separate domains in the complex plane.Bifurcation of the mode eigenvalue (and the ensuing mode conversion) occurs when the toroidal rotation frequency and thermal transport coefficients exceed certain critical values-in our case a rotation frequency between ω 0 /ω A = 0.0056 and 0.0057, and the thermal transport coefficient between χ ⊥ = 8 × 10 −5 and 9 ×10 −5 (χ ∥ /χ ⊥ = 10 6 is fixed).The behaviour of the two IK branches is qualitatively different in the complex plane.For instance, at χ ⊥ = 10 −4 (labelled by yellow circles), two branches are either stabilized or destabilized with increasing plasma flow speed.This is the same as those shown in figure 3(a).However, at fixed χ ⊥ = 8 × 10 −5 (magenta triangles), one branch is first stabilized by increasing plasma flow, then followed by destabilization, while the other branch experiences the opposite trend as the plasma flow speed increases from ω 0 /ω A = 0.005-0.006.
We point out that the aforementioned mode conversion also depends on the ratio of the parallel to perpendicular transport coefficients χ ∥ /χ ⊥ .This ratio has so far been fixed at 10 6 .Figure 6 reports results with even higher values for this ratio.As χ ∥ /χ ⊥ increases (at fixed plasma rotation), the branch-1 mode is strongly destabilized while branch-2 is weakly stabilized.This opposite effects on the mode growth rate tend to separate the two branches, leading to disappearance of the mode conversion phenomenon.
In conclusion, we have reported that anisotropic thermal transport, together with plasma flow, can fully stabilize the IK but via two completely different mechanisms.One is the synergetic stabilization due to high level of anisotropic thermal transport and plasma flow.Full stabilization in this case requires fast (but still subsonic) plasma flow, with a threshold value of ω 0 /ω A = 0.0085 when χ ⊥ = 10 −4 in our example.The other stabilization mechanism involves indirect effect of thermal transport, with the latter induces finite mode frequency in the plasma frame thus allowing the mode to be in resonance with the stable continuum shear Alfvén waves.A complete stabilization of the mode can also be achieved in this case even in a static plasma.This kind of stabilization depends on the resonance matching condition as shown in equation ( 7), which is not always possible to be satisfied.For instance, in the absence of anisotropic thermal transport, the IK generally rotates with the plasma and the mode frequency in the plasma frame is too slow to be in resonant with plasma continuum waves.This is also why such a stabilization mechanism has not been previously reported.Inclusion of anisotropic thermal transport (and maybe other new physics effects beyond the scope of the present work) makes this kind of continuum wave resonance possible, and consequently the access to mode stabilization.Moreover, the combined effects of plasma rotation and thermal transport introduce the possibility of mode eigenvalue bifurcation and mode conversion (at relatively low levels of thermal transport anisotropy) for the IK.The findings reported in this work thus provide new physics insights in the IK stability in tokamak fusion plasmas.
Although not the focus of the present work, we point out that non-linear effects associated with IK, typically leading to sawtooth, also play important roles in the plasma dynamics.Non-linear MHD simulations, such as that reported in [26], are required to quantify the sawtooth behaviour, which is outside the scope of the present study.Furthermore, physics beyond the MHD model, such as kinetic effects associated with fast ions, are known to be also important for the sawtooth dynamics [6,27].In particular, fusion-born α-particles can strongly stabilize the IK, resulting in monster sawtooth with a potential risk of triggering large neoclassical islands and plasma disruption in future tokamak reactors.The destabilizing effect of anisotropic thermal transport on the IK at high values of χ ∥ /χ ⊥ , that we identify in this study, may be favourable for avoiding monster sawtooth by counteracting the fast ion stabilization.

Figure 1 .
Figure 1.Toroidal equilibria, from the HL-3 tokamak, assumed in MARS-F modelling showing (a) the plasma boundary and wall shapes (b) two equilibrium safety factor profiles with the on-axis values of q 0 = 0.90 (dashed) and q 0 = 0.97 (solid), (c) two equilibrium pressure profiles normalized by B 2 0 /µ 0 , and (d) the toroidal rotation frequency normalized by the toroidal Alfvén frequency ω A .The radial coordinate s ≡ ψ 1/2 p , with ψ p being the normalized equilibrium poloidal magnetic flux.

Figure 2 .
Figure 2. The MARS-F computed growth rate of the n = 1 internal kink mode with varying on-axis safety factor, assuming (a) vanishing, (b) a finite (ω 0 /ω A = 0.01) toroidal rotation of the plasma.Compared are cases without the thermal transport terms χ ⊥ = 0 (w/o TTT, blue solid lines), with finite thermal transport χ ⊥ = 10 −5 (black dashed lines) and χ ⊥ = 10 −4 (red dashed line), respectively.The ratio of the parallel to perpendicular thermal transport coefficients, χ ∥ /χ ⊥ , is kept fixed at 10 6 for cases with finite thermal conduction.The value of χ ⊥ is normalized by a 2 /τ A , with a being the plasma minor radius and τ A being the toroidal Alfvén time.

Figure 3 .
Figure 3.The computed (a) growth rate and (b) real frequency of the n = 1 internal kink with varying toroidal rotation frequency while fixing the thermal transport coefficients (χ ⊥ = 10 −4 , χ ∥ /χ ⊥ = 10 6) and the on-axis safety factor ( q 0 = 0.97).Shown are two branches of the mode, indicated by red-solid and blue-dashed lines, respectively.

Figure 4 .
Figure 4.The m/n = 1/1 component of the radial displacement associated with the internal kink mode, assuming (a) a finite (ω 0 /ω A = 0.015 corresponding to branch 1 shown in figure3) and (b) vanishing plasma flow (corresponding to branch 2 from figure3).The latter is also plotted along the safety factor q in the lower panel of (c), together with the shear Alfvén wave frequency k ∥ v A (upper panel) near the q = 1 surface (indicated by vertical dashed lines in all plots).The horizontal line in the upper panel of (c) indicates the computed mode frequency in the plasma frame.The two solid vertical lines in (c) indicate the radial locations where the mode frequency matches the shear Alfvén wave frequency.Assumed is a case with finite thermal transport coefficients (χ ⊥ = 10 −4 , χ ∥ /χ ⊥ = 10 6 ) and an equilibrium with the on-axis safety factor of q 0 = 0.97.