Simulations of nonlinear interaction between beta-induced Alfvén eigenmode and tearing mode

The nonlinear interaction between the m/n=2/1 beta-induced Alfvén eigenmode (BAE) and the m/n=2/1 tearing mode (TM) observed in the HL-2A experiment is systematically studied with the hybrid kinetic-magnetohydrodynamic code M3D-K. It is found that the nonlinear growth of TM can lead to a gradient buildup on the magnetic island (MI) edge, and then triggers the destabilization of the linearly-stable BAE. Then, the triggered BAE together with TM can produce a significant redistribution of energetic particles. For BAE linearly-dominant-unstable case, the TM activity results in the inward motion of BAE mode structure, and the BAE can have a delay effect on the MI saturation. Furthermore, a high-frequency axisymmetric m/n=0/0 ‘breathing’ mode is generated by the mode coupling of BAE and TM, agreeing well with the experimental observation, and causes the synchronized periodic oscillation of MI width.


Introduction
Tearing mode (TM) [1,2], one of the most dangerous magnetohydrodynamic (MHD) instabilities in tokamak plasmas, can lead to the formation of magnetic islands (MIs), increase local radial transport and degrade plasma confinement, even trigger disruptions. The TM can also interact * Author to whom any correspondence should be addressed.
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. with energetic particles (EPs) and leads to the EP redistribution/losses [3][4][5][6]. Hence, the suppression or control of TM is critical to plasma confinement and has been intensively studied [7][8][9][10]. On the other hand, the Alfvén eigenmodes (AEs), destabilized by EPs and extensively observed in magnetic fusion devices, are crucial to the performance of burning plasmas [11][12][13]. For example, the beta-induced AE (BAE) [14,15], located in the gap caused by the combination of finite plasma beta and geodesic curvature, can be driven by EPs through wave-particle resonance, and can cause the EP redistribution/losses in turn. A simulation of the China Fusion Engineering Test Reactor (CFETR) hybrid scenario suggests that the AEs can reduce the EP density in the core and induce approximately 30% of EP transport from the core region to the edge [16]. Moreover, the EP redistribution/losses not only reduce the heating efficiency and deteriorate the fusion reactor performance, but also contribute significantly to plasma-wall-interaction due to the large power carried by them [17,18]. Based on the interaction of EPs with a single AE, further investigations have advanced to the generation of AE spectra majorly due to mode-mode coupling, particularly in the burning plasma regime, and their backreaction on transport coefficients of EPs [11,19].
Experiments in DIII-D tokamak have suggested that the synergistic effect of AEs and TM is likely responsible for the intensified EP redistribution/losses [20]. Moreover, the nonlinear interactions between AEs and TM, especially for the nonlinear interaction of BAE and TM, are often observed in many tokamak experiments [21][22][23][24][25][26][27][28]. Experiments in JET and EAST have reported that the nonlinear interaction between TM and geodesic acoustic mode is responsible for the excitation of twin BAEs [24,26]. Besides driving BAE, EPs can also suppress or destabilize TM. TM and BAE can modulate each other, and cause the redistribution/losses of EPs [29]. Thus, theoretical investigation [30][31][32] and numerical simulations [4,33] are also carried out to explore the profound physical mechanism. Clearly, the interaction between AEs, TM and EPs is so complex that shows comprehensive approaches to wave-particle interaction and EPs confinement improvement.
Since the mode structures of BAE and TM are both localized around the rational surfaces q = m/n (m and n are the poloidal and toroidal mode numbers), they can be expected to interact with each other strongly. In a recent HL-2A experimental investigation [27], it has been indeed found that with prominent TM activity, the nonlinear interaction of BAE and TM generates a high-frequency axisymmetric n = 0 mode by three-mode coupling and causes significant redistribution/losses of EPs.
Based on the experiment in HL-2A, we then perform a numerical simulation to explore the underlying mechanism of the nonlinear interaction referring to BAE, TM and EPs. We have reproduced the experimental observation that a highfrequency axisymmetric m/n = 0/0 'breathing' mode (BM) is generated by the mode coupling of BAE and TM. The interaction process shown in figure 1 is revealed and described as follows. By magnetic free-energy stored in the bulk plasma current, a TM is generated and observed. Besides heating/ accelerating the bulk plasma as well known, the growth of TM-MI also significantly changes the EPs distribution to trigger BAE. Then, the triggered BAE together with TM greatly enhance the redistribution of EPs. The BM for completing the three-mode coupling cycle also modulates the MI and leads to a synchronized periodic oscillation of m/n = 2/1 MI width, different from that caused by the sawtooth crashes [10,34] and the plasma rotation [35]. These results are novel and could be important for current tokamak experiments and future burning plasmas such as ITER and CFETR.
This paper is organized as follows. The M3D-K code and simulation parameters used in this work are briefly described in section 2. The simulation results are presented in section 3. Finally, the discussion and conclusion are given in section 4.

Experimental phenomena and simulation setup
The experimental result performed in deuterium plasma is presented in figure 2, where MHD activities in HL-2A shot #30759 are shown in panel (a). It can be found that the frequencies satisfy the relation: ω BM = ω BAE − ω TM , revealing a typical three-mode coupling process of BAE, TM and BM. The corresponding safety factor and total/EP β profiles are plotted in panels (b) and (c), where β is the ratio of plasma pressure to magnetic pressure. The global kinetic-MHD hybrid code M3D-K is employed in this work. Here, the hybrid model is a combination of MHD model and kinetic model, and it is more powerful than MHD model and more economical than full-kinetic model, although the kinetic effects of bulk plasma, such as Landau damping and finite Larmor radius, are not included. In M3D-K, the bulk plasma is described by the resistive MHD equations and the energetic ions (EPs hereafter) are described by the drift-kinetic equations, and the kinetic effects of EPs are coupled into the momentum equation through the EPs pressure tensor P h [36,37]. This hybrid model is widely applied to study the problems of the interaction between EP and MHD instabilities [3][4][5][38][39][40][41][42][43][44][45]. The main parameters are as follows. Major radius R 0 = 1.65 m, minor radius a = 0.4 m, toroidal magnetic field B 0 = 1.35 T, central number density n 0 = 2.0187 × 10 19 m −3 , Alfvén speed v A = B 0 /(µ 0 ρ 0 ) 1/2 = 4.63 × 10 6 m s −1 and Alfvén time τ A = R 0 /v A = 3.56 × 10 −7 s, central total beta β total = 0.9% including both thermal plasma β thermal and EPs β h , the EPs' pressure fraction β frac = β h /β total = 0.50. The injection energy of NBI is E 0 = 45 keV, and the corresponding injection velocity of where c is a normalization factor, which is chosen to match the peak beta, H is the step function, v 0 is the cutoff velocity of EPs, and v c = 0.58v 0 is the critical velocity. Λ ≡ µB 0 /E is the pitch angle, Λ 0 = 0.0 is the central pitch angle, ∆Λ = 0.30 is the width of pitch angle distribution. ∆ψ is radial width for the distribution, ⟨ψ⟩ is ψ averaged over the EPs orbit. Based on the experimental observation and for simplicity, we only retained the toroidal mode numbers of n = 0, 1 in simulation to investigate the interaction of m/n = 2/1 BAE and m/n = 2/1 TM more clearly, and all simulations are carried out with both n = 0 and n = 1 modes.

Pure simulations of m/n = 2/1 BAE and m/n = 2/1 TM
Firstly, the simulation results of TM-only cases (without EPs) and BAE-only cases (resistivity η = 0) are performed in this subsection. For the given equilibrium parameters, the TM is stable without resistivity, and the BAE is stable without EPs. For TM-only cases, as shown in figure 3(a), the scaling law of linear growth rate and resistivity satisfies as γ ∝ η 3/5 , which is consistent with the previous result and analytical theory [33,46]. For BAE-only cases, from the green curve in figure 3(b), it can be clearly found that widening EP distribution by raising ∆ψ reduces the gradient of the distribution and stabilizes BAE. The linear eigenmode structures of stream function (U, related to the incompressible part of the plasma velocity) are plotted in 3(c) and (d), one can see that both TM and BAE exhibit a clear m/n = 2/1 mode structure. The TM is localized around the q = 2 surface, while the BAE is mainly located inside the q = 2 surface.

Influence of TM activity on BAE
In this subsection, the BAE + TM simulations are performed with resistivity η = 4 × 10 −5 for various EP distribution radial width ∆ψ. As shown in figure 4, the linear growth rate of TM slowly rises with the increased EPs β h at q = 2 surface, consistent with the previous simulation [3]. When ∆ψ = 0.20, the EP-driven BAE is the linearly-dominant-unstable mode because its growth rate is much larger than that of the TM. The U cos,n1 is decomposed by a high-pass filter into two contributions of the low-frequency TM (ω ⩽ 0.05ω A , red line) and the high-frequency BAE (ω > 0.05ω A , blue line). The corresponding frequency spectrum is plotted in figure 5(c). One can see that the BAE component is triggered after the TM saturation around t = 1500τ A , and finally saturates from t = 2500τ A with a weakly downward frequency chirping.
The mode-particle interaction in the process is presented by the resonance condition of nω ϕ + pω θ = ω, where n is the toroidal mode number, p is an integer. ω, ω ϕ and ω θ are the mode frequency, the toroidal and poloidal transit frequency for passing ions. Figure 6 shows the Alfvén continuous spectrum calculated by NOVA code [47] and the perturbed distribution function (δf). Firstly, as shown in figures 6(a) and (d), only the low-frequency TM exists at t = 1500τ A , and the broad hole-clump structure is related to the p = −2 TM resonance, corresponding to the redistribution of EPs induced by the 2/1 MI. Then, at t = 2000τ A , the high-frequency BAE emerges with frequency being located inside the BAE-gap in Alfvén continuum (see figure 6(b)), agreeing well with  experimental measurements of f BAE = 0.125ω A /2π + f rot ∼ 67.9 kHz, where f rot is the plasma rotation frequency. As shown in figure 6(e), below the TM-induced hole-clump pair, a new hole-clump pair corresponding to the p = −1 BAE resonance appears simultaneously with the destabilization of BAE. Finally, as shown in figure 6(c) and ( f ), the lowfrequency TM and the high-frequency BAE coexist in the late nonlinear stage. The hole-clump pair of p = −1 BAE resonance becomes significant, and the resonance regions move to the lower energy E with the clump structure being dominant, corresponding to the weakly downward frequency chirping of BAE [38].
To further elucidate the excitation mechanism of BAE by TM, the radial profiles of distribution f(P ϕ ) and safety factor q at different time are plotted in figure 7(a). The analysis results can be shown as follows. (a) The q profile is flattened around the q = 2 surface after the saturation of TM, and so is f(P ϕ ) due to EP transport induced by the MI. (b) The gradient df/dP ϕ near the 2/1 MI separatrix increases. (c) Meanwhile, the core β h gets larger than the edge β h , due to the spatial exponential distribution chosen in equation (1). Thus, the combination of factors (a)-(c) destabilizes the linearlystable BAE near the 2/1 MI separatrix and inside the initial q = 2 surface in the nonlinear stage. Correspondingly, a BAE-EP interaction induced flattening region of f(P ϕ ) emerges in the late stage after t = 2000τ A . For BAE linearly-dominantunstable case (∆ψ = 0.20), the above three factors can lead to the inward motion of BAE structure, as shown in figure 8. Moreover, the EP redistribution level g is measured by the deviation of the distribution f (t) from its initial f(t = 0), as   Here the toroidal angular momentum P ϕ can be regarded as a radial variable. The maximal P ϕ of EPs around E = 0.51E 0 is about −0.346, and the EPs with P ϕ > −0.346 are lost by the prompt loss [56]. The red region denotes the radial location of BAE mode structure, and the red curves represent the q profiles at t = 0τ A and t = 1500τ A . (b) Time evolution of redistribution level g of EPs.  g =´| f(t) − f(t = 0)|dP ϕ /´f(t = 0)dP ϕ , and its time evolution is plotted in figure 7(b). One can see that in the early stage, the level g is mostly affected by the radial profile flattening due to the topology change inside the growing 2/1 MI, and then slowed down due to the TM saturation. Then, the redistribution level increases significantly after the excitation of the BAE to a level where the BAE together with the TM lead to the enhanced EP redistribution in the late nonlinear stage.

Mode coupling between BAE and TM
To reveal the three-mode interaction process, the evolution of U cos,n0 and its corresponding frequency spectrum are plotted in figures 9(b) and (c). The amplitude of U cos,n0 exhibits a second growth and a high-frequency oscillation after the excitation of BAE around t = 1500τ A , where the former corresponds to the generation of the n = 0 component of BAE, and the latter indicates the excitation of the high-frequency axisymmetric m/n = 0/0 BM. Figure 9(c) shows the evolution of the perturbed parallel current ∆C BM for the BM, with ∆C = −R∆J ϕ , ∆J ϕ is the perturbed toroidal current. The axisymmetric BM at midplane (Z = 0) exhibits a periodic oscillation during the generation of three-mode coupling (after t = 1500τ A ) with a period T = 2π/ω BM ∼ 70τ A ∼ 25 µs. For the three coupled modes of BAE, TM and BM in the nonlinear stage, both matching conditions of the mode number and the frequency are satisfied [48] as (n, m) BM = (n, m) BAE − (n, m) TM and ω BM = ω BAE − ω TM (see red curve in figure 9(b)), consistent with the experimental observation. Therefore, this BM is generated in the nonlinear threemode coupling by the high-frequency BAE and the lowfrequency TM.

Effects of BM and BAE on 2/1 MI evolution
In this subsection, the effects of BM and BAE on 2/1 MI evolution are investigated respectively. As shown in figure 10(a), the saturation width of 2/1 MI is about 0.2a ∼ 8 cm, and it is intriguing to find that the m/n = 2/1 MI width exhibits a periodic oscillation after the generation of BM with a period of T ∼ 25 µs, the same as that of the BM. Clearly, it is a BM modulation on TM. It is found that this oscillatory MI width has a variation rate of ∼10%. In comparison with the 2/1 MI width, the 3/1 MI width is much smaller and remains nearly unchanged after its saturation. Furthermore, the BM toroidal current ∆C BM driven by the mode coupling of 2/1 BAE and 2/1 TM is localized around the 2/1 MI separatrix, and changes periodically as shown in figure 10(b) and (c), corresponding to a periodic current modification.
Then, the mechanism of this MI width oscillation is further investigated. Figures 11(a)-(c) shows the evolution of three high-frequency branches of BAE, BM and W HF , with  W HF being the high-frequency oscillation component of the MI width by filtering out the low-frequency (ω ⩽ 0.05ω A ) contribution to the MI. Their relation is then measured by the Pearson correlation coefficient method [49,50] P(x, y) = where P(x, y) describes the linear correlation between the variables x and y with −1 ⩽ P ⩽ 1. The larger the magnitude of P(x, y) is, the stronger the correlation between x and y is.  [45] shown in figure 11(d), and their phases are locked at 7π/4. From the enlarged view of the temporal evolution of ∆C BM as illustrated in figure 9(c) one then again clearly sees that the highfrequency oscillation of the MI width is synchronized with the periodic change of ∆C BM . In tokamak plasmas, the evolution of MI depends on the current profile [51]. One can also say that the nonlinear oscillation of 2/1 MI width is induced by the periodic current modification of ∆C BM localized around the 2/1 MI separatrix. Finally, the effect of BAE on 2/1 MI evolution is investigated in figure 12. In comparison with the TM-only case, the MI saturation time remains almost unchanged for TM linearly-dominant-unstable case (η = 9 × 10 −5 ) as well as BAE and TM linearly-coexisted case (η = 4 × 10 −5 ). However, the MI saturation time is prolonged for BAE linearlydominant-unstable case (η = 1 × 10 −5 ), which is similar to the MI saturation time delay induced by the external driven current [52,53]. The delayed saturation time of MI by BAE is beneficial to control TM, which may be associated with the current modulation of the n = 0 component of BAE, and it will be further studied in future work.

Discussion and conclusion
Based on the experimental observation in HL-2A, we performed a systematic study of the nonlinear interaction between m/n = 2/1 BAE and m/n = 2/1 TM by using the kinetic-MHD code M3D-K. In the presence of TM, the linearly-stable BAE can be triggered due to the combined effects of flattened q profile at q = 2 surface and steepened EP radial gradient caused by the MI-induced EP redistribution in the nonlinear stage. This combination can also lead to the inward motion of BAE structure for BAE linearly-dominant-unstable case. The BAE together with the TM can lead to the excessive EP redistribution. It is also found that the BAE can have a delay effect on MI saturation, which is beneficial for MI control, and future work is needed to fully unravel its mechanism.
Moreover, we have reproduced the experimental observed high-frequency axisymmetric m/n = 0/0 BM, which is generated by the nonlinear mode coupling of BAE and TM. The perturbed current of BM (∆C BM ) is localized around the 2/1 MI separatrix, with a periodic amplitude oscillation in a period T ∼ 25 µs. Then, the periodic current modification from the ∆C BM produces the synchronized periodic oscillation of 2/1 MI width, with a nonlinear disturbance rate of ∼10%. The 3/1 MI width remains almost invariant without BM around q = 3 surface in the nonlinear stage.
Besides affecting the plasma transport [54] as well as regulating the turbulence [55] and current in and out of MI, the high-frequency axisymmetric n = 0 BM can also provide an energy channel between different space-time scales [27]. Furthermore, the periodic oscillation of MI width caused by a high-frequency axisymmetric n = 0 BM can be ubiquitous in fusion plasma. The oscillation of MI width could induce neighboring MIs interaction, leading to the stochasticity of magnetic field to trigger the plasma disruption. These new findings allow us to better understand the underlying interaction physics and develop relevant predictions to prevent EPs redistribution/losses and improve plasma performance.