Waveform distortion of off-axis fishbone instability in the nonlinear magnetohydrodynamic evolution in tokamak plasmas

We have investigated the waveform distortion of energetic particle driven off-axis fishbone mode (OFM) in tokamak plasmas with kinetic magnetohydrodynamic (MHD) hybrid simulations. We extended our previous simulations (Li et al 2022 Nucl. Fusion 62 026013) by considering higher-n harmonics in the MHD fluid, where n is toroidal mode number. The waveform distortion is successfully reproduced in the simulation for both magnetic fluctuations and temperature fluctuations. It is clarified that the waveform distortion arises from the superposition of the n = 2 harmonics on the fundamental n = 1 harmonics of OFM, where the n = 2 harmonics are generated by the MHD nonlinearity from the n = 1 OFM. Two types of waveform distortion can occur depending on the phase relationship between the n = 1 and n = 2 harmonics and the relative amplitude of the n = 2 harmonics to the n = 1 harmonics. Lissajous curve analyses indicate that the wave couplings between the n = 1 and n = 2 harmonics with phase-lock ∼π and ∼0 lead to ‘rising distortion’ and ‘falling distortion’, respectively. The two types of waveform distortion can be attributed to the strong shearing profile of radial MHD velocity with n = 2 around the q = 2 magnetic flux surface. The dependence of waveform distortion on viscosity is investigated. It is found that the viscosity which is needed to reproduce the waveform distortion is larger than that in the experiment.


Introduction
Off-axis fishbone mode (OFM) is an energetic-particle (EP) driven magnetohydrodynamic (MHD) instability, also called EP driven wall mode. The OFMs have been observed in tokamak plasmas such as JT-60U [1], DIII-D [2,3], JET [4] and NSTX [5] plasmas with the central safety factor (q) well above unity. Trapped EPs destabilize OFMs, and the frequency chirps down after the saturation of the mode amplitude growth. These properties are similar to those of classical q = 1 fishbones [6]. OFMs, on the other hand, are located close to the q = 2 magnetic flux surface, whereas classical fishbones are inside the q = 1 surface. OFMs are essential for fusion research because they redistribute EPs, leading to the destabilization of resistive wall modes (RWMs) [7] and edge localized modes (ELMs) [8].
The mechanisms underlying interactions between OFM and other MHD modes with marginal stability have been studied theoretically and experimentally. During OFM bursts in DIII-D, the substantial reduction of both the density of energetic particles [9] and plasma rotation frequency [10] are supposed to contribute to the destabilization of the RWM. Since the instantaneous change in neutron emission and the toroidal velocity on the q = 2 surface was confirmed, some researchers proposed that nonambipolar radial electric field accumulation could be a possible mechanism for rotation decreasing [11,12]. A theoretical work revealed that a mode conversion could occur between the RWM and precessional OFM [13]. It was also reported that the OFM interacts with ELMs on JT-60U [14,15], resulting in the ELMs being triggered by the OFM in high-β operation. It was inferred that the additional pressure caused by the EP transport by the OFM changes the edge stability and triggers the ELM.
In a previous work, we conducted a self-consistent simulation study of OFM [16]. We showed the spatial structure of OFM with a strong shearing shape and the resonance condition with deeply trapped energetic particles. In the nonlinear stage, we found that the resonance frequency of the particles that transfer energy to the OFM chirps down with the split of resonant regions in phase space, which may result in the chirping down of the OFM frequency. At the same time, the distribution function is flattened along the E ′ = const. line, leading to the saturation of the instability. E ′ is a combination of energy and toroidal canonical momentum, and conserved during the wave-particle interaction neglecting the time variation of the mode amplitude and frequency [17][18][19].
The most exciting feature of OFM in experiments is the nonlinear mode distortion at the maximum mode amplitude. The amplitude of OFM increases exponentially in a nearly sinusoidal waveform. At the peak amplitude, the waveform deviates significantly from a sine wave [12,14]. Measurements with the DIII-D and JT-60U toroidal arrays show that the increase in mode distortion occurred concurrently with the increase in higher-n components and energetic particle loss, where n is toroidal mode number. One possible mode distortion in high-beta plasmas was reported by Park [20] with n = 1 distortion in TFTR. Once the high-beta plasma reaches the ideal mode limit, a low n mode is destabilized initially, leading to nonlinear evolution of the plasma in 3D equilibrium. The resulting equilibrium exhibits a steep pressure gradient in the bad curvature region. This localized pressure gradient can then induce a new high-n mode, causing a local bulge in pressure. In contrast to the OFM case, where the mode distortion is spread throughout the mode structure, the distortion in the TFTR case is localized on the steep pressure gradient. Therefore, it is unlikely that the ballooning effect contributes to the OFM distortion.
In this paper, we will present a new mechanism where MHD nonlinearity causes the waveform distortion of OFM. We employ hybrid simulations of MHD fluid interacting with energetic particles. The simulations reproduce the waveform distortion of OFM for the first time, which is primarily caused by MHD nonlinearity generating n = 2 harmonics. Thus, the amplitude and the spatial structure of n = 1 OFM play essential roles in the formation of distortion. In this paper, we will show that there are two types of waveform distortion, rising distortion and falling distortion, depending on the phase relationship between the n = 1 OFM and n = 2 harmonics as well as the relative amplitude of the n = 2 harmonics to the n = 1 OFM. In addition, the Lissajous figures show the phase-lock between n = 1 OFM and n = 2 harmonics for these two types of waveform distortion. Through the analysis of nonlinear components in MHD equations, the difference between the two types of waveform distortion is attributed to the shearing structure of the n = 2 component of radial MHD velocity v rad,n=2 . This paper is organized into the following sections. In section 2, the simulation model is described. In section 3, the simulation results are presented, focusing on the nonlinear behaviors of OFM. Section 4 is devoted to discussion and summary.

Simulation model
The same code for kinetic-MHD hybrid simulations, MEGA [21][22][23][24], was used as the previous work [16]. We use the whole tokamak plasma domain with the toroidal angle range 0 ⩽ ϕ < 2π. Then, the simulation region is R c − a ⩽ R ⩽ R c + a, 0 ⩽ ϕ < 2π and −1.7a ⩽ Z ⩽ 1.7a, where R c = 1.7 m and a = 0.6 m are the major radius and the minor radius, respectively. The number of grid points is 128 × 64 × 128 for cylindrical coordinates (R, ϕ, Z) with approximately 8 × 10 6 particles. When analyzing the simulation results, magnetic flux coordinates (ψ, θ, ϕ) are used, where ψ is flux surface label, ϕ is toroidal angle, and θ is poloidal angle. The values of viscosity and diffusivity are set to be ν = 10 −6 v A R c and ν n = χ = 0, and the resistivity is η = 10 −7 µ 0 v A R c in the simulations unless otherwise specified, where v A is the Alfvén velocity at the plasma center.
An anisotropic slowing-down distribution is used here to be consistent with the distribution of energetic particles in the experiments. The distribution function is given by whereψ is normalized poloidal magnetic flux withψ = 0 at the plasma center andψ = 1 at the plasma edge, and ∆ψ = 0.3. Pitch angle variable is represented by Λ = µB 0 /E, where µ, B 0 , and E are magnetic moment, the magnetic field strength at the plasma center, kinetic energy of particles with Λ 0 = 1.1 and ∆Λ = 0.1 being the distribution peak location and width, respectively. The background plasma is deuterium with the number density of 3 × 10 19 m −3 . The magnetic field strength at the plasma center is B 0 = 1.7 T. The injection velocity of energetic particles is v inj = 0.58v A , corresponding to 80 keV deuterium neutral beam with ∆v = 0.1v A . The critical velocity is v crit = 0.62v A . The profiles of EP beta, bulk plasma beta and safety factor are shown in figure 1. The bulk plasma and EP beta profiles are defined by with β 0 the value at the center. Regarding the safety factor, the on-axis value and the edge value are q r=0 = 1.6 and q r=a = 5.0.

The basic mechanism of waveform distortion
Off-axis fishbone bursts generally last 2-3 ms and occur approximately every 20 ms in the experiments [1,3,12,25]. A typical Mirnov coil signal has a nearly sinusoidal waveform in the initial phase, and the amplitude grows roughly exponentially. Around the mode saturation, the waveform departs significantly from a sinusoidal wave. Simultaneously, the toroidal array detects the emergence of higher-n harmonics. This novel nonlinear property was reproduced for the first time in our simulations. We extend the simulation retaining high-n MHD harmonics with n > 1 while they were filtered out in our previous work. The waveform distortion for both magnetic and temperature fluctuations is successfully observed in our simulations. Additionally, some cases involving higher-n harmonics up to n = 7 have also been studied, as shown in figure 2, and it is clear that the n = 2 component is the primary contributor to the occurrence of waveform distortion. The time evolutions of the MHD perturbation energy (W) with n = 1, 2 and the frequency of radial MHD velocity with m/n = 2/1, 4/2 are shown in figure 3. Both the n = 1 and n = 2 harmonics have exponential growth and the saturation of the MHD perturbation energy at t = 1.68 ms. The growth rates of radial MHD velocity are 0.0046 ω A and 0.0092 ω A for n = 1 and n = 2 harmonics, where ω A is the Alfvén frequency at the plasma center, while the mode frequencies are   ω 2/1 = 11.2 kHz and ω 4/2 = 22.1 kHz before the saturation and they both chirp down approximately 40% along with the mode damping. The n = 1 OFM is driven by deeply trapped particles, with the mode frequency being equal to the precessional drift frequency of energetic particles. The growth rate and mode frequency of the n = 2 harmonics are twice those of the n = 1 harmonics. We can conclude that the n = 1 harmonics generate the n = 2 harmonics through the MHD nonlinearity. Next, we will explain how the n = 2 harmonics contribute to the formation of waveform distortion.
In the simulation, δṪ e = dδT e /dt is investigated first to compare with the δT e −ECE signal in the experiment. Figure 4 shows the wave oscillation measured on the q = 2 magnetic surface in the outer midplane and poloidal structures for n = 1, n = 2 and total component. In figures 4(i), the oscillations of n = 1, 2 harmonics maintain the sinusoidal structure, while the total component exhibits experiment-like distortion. This clearly indicates that the superposition of the n = 2 harmonics on the n = 1 harmonics of OFM causes the waveform distortion. The observed OFM in experiments rotates evenly in the toroidal direction, which is consistent with our simulation results shown in figure 4(i). In figures 4(iii) and (iv), snapshots of the spatial structure of δṪ e on the poloidal plane with ϕ = 0 are given for both the linear and the nonlinear phases. Before the saturation at t = 1.23 ms, the n = 1 harmonics are composed of m = 2 structure around the q = 2 magnetic surface, as indicated by the bold solid line. For the n = 2 harmonics, they have m = 4 outside the q = 2 surface, while m = 3 is dominant inside the q = 2 surface. The generation of the m/n = 3/2 harmonics can be attributed to the nonlinear interaction between the 1/1 and 2/1 harmonics. Changing the bulk plasma profile from uniform [16] to Gaussian causes an increase in the amplitudes of the 1/1 and 2/1 harmonics, which in turn leads to a strong nonlinear generation of 3/2 harmonics. In particular, we notice the shearing characteristic of the n = 1, 2 harmonics on the poloidal plane. The nonperturbative kinetic effect of energetic particles is one probable mechanism for the shearing structure of the OFM [16]. After the saturation, the n = 2 mode grows to a comparable amplitude to the n = 1 mode, resulting in a stronger shearing profile of the total component as illustrated in figures 4((c)-(iv)), which can also be seen in figures 10((b)-(iii)) and ((d)-(iii)). Based on these analyses, we can confirm that the n = 2 harmonics substantially contribute to the waveform distortion.
The coupling effects between the n = 1 and n = 2 harmonics are examined in depth here to understand the mechanism of waveform distortion better. The temporal evolution of radial profiles of δṪ e for n = 1, n = 2 and total components are presented in figures 5(a)-(c), respectively. Furthermore, the waveforms of δṪ e measured at r/a = 0.35 and 0.64 are plotted on the figures. In figure 5(a), the oscillation is located at 0.1 ⩽ r/a ⩽ 0.67 after the saturation at t = 1.68 ms. The n = 2 mode generated by MHD nonlinearity oscillates primarily in the range of 0.17 ⩽ r/a ⩽ 0.4 and 0.4 ⩽ r/a ⩽ 0.65 as shown in figure 5(b), while the shearing structure is located at r/a ∼ 0.4. We should emphasize that these two regions have a phase difference of nearly π caused by mode shearing. The phase difference is a critical factor in waveform distortion formation, as we will discuss later. It is found that the waveforms at r/a = 0.35 and 0.64 for n = 1 and n = 2 remain sinusoidal even in the late stage of simulation. Nonetheless, as illustrated in figure 5(c), the total δṪ e component has strong coupling effects, with n = 2 harmonics significantly impacting the radial structure of OFM. Consequently, we have observed two types of waveform distortion located around r/a = 0.35 and 0.64. Based on the location of distortion in the waveform, we categorize these into 'falling distortion' and 'rising distortion' as shown in figures 5(d) and (e). In particular, the rising distortion is formed by two sinusoidal waves with a phase difference of π and is found at r/a ≃ 0.6, which is identical to the experimentally observed waveform distortion [3,12,14].
The type of waveform distortion depends on the phase relationship. The Lissajous figure can help us understand the phase relationship between the n = 1 and n = 2 harmonics since the specific curves indicate the locked phase between two waves [26,27]. Figures 6(a) and (c) show δṪ e and the envelope of δṪ env e for n = 1 and n = 2 harmonics at r/a = 0.35 and 0.60, and the signals are normalized to δṪ e /δṪ env e to extract the phase as shown in figures 6(b) and (d). The phase relationship between the n = 1 and n = 2 harmonics is locked at ∼−10 (∼ 170) deg at r/a = 0.35 (r/a = 0.64). Theoretical and numerical studies have shown that phase-locking could result in meso-or macro-scale secular EP transport [28]. Notably, the locking phase has a good agreement with the diagram shown in figures 5(d) and (e).

Nonlinear components analysis
The nonlinear terms in the MHD equations generate fluctuations with toroidal mode number of n = 2 in the evolution of temperature fluctuations of OFM. It would be interesting to learn which factor is critical for the development of n = 2 harmonics, and may result in the two types of waveform distortion at different locations. The temperature evolution is given by The dissipation terms can be neglected in the linear phase, when we analyze the sideband harmonics [23]. Thus, the δṪ e,n=2 primarily consists of −v n=2 · ∇T eq , −v n=1 · ∇T n=1 and −(γ − 1)T eq ∇ · v n=2 as shown in figure 7. It is found that the −v n=2 · ∇T eq term is the dominant component, in which the spatial profile of v n=2 is crucial to the evolution of the n = 2 harmonics of temperature fluctuations. Because the equilibrium temperature profile is uniform in the ϕ and θ directions, only v rad,n=2 has a significant effect on the spatial profile of δṪ e . We see in figure 8 that the primary component of v rad,n=2 is m = 4 and the phase which is the ratio of the sine component to the cosine component changes abruptly at r/a ∼ 0.4. The phase of v rad,n=2 distinguishes falling and rising distortion.

Effects of weak dissipation
Since the applied MHD model lacks kinetic damping, we employ artificial dissipation coefficients to control the damping. The viscosity and resistivity coefficients employed in this work are ν = 10 −6 v A R c and η = 10 −7 µ 0 v A R c , as we mentioned in section 2. Typical dissipation coefficients in real fusion plasmas are lower by orders of magnitude than simulations. Thus, it is important to examine the effects of weak dissipation. We performed several runs for weak viscosity and resistivity, and found that the mode structure and the saturation level are insensitive to the resistivity coefficient, while the viscosity coefficient plays an important role. The saturation level of mode energy of n = 1 harmonics decreases by 10% in the weak viscosity simulation, while in contrast, the energy of n = 2 harmonics has been doubled. The variation in the saturation amplitude is attributed to the nonlinear MHD effects with weaker dissipation [23,24]. In addition, the experimental waveform distortion of OFM is difficult to be observed for lower viscosity value as shown in figure 9, where the n = 2 harmonics has a larger amplitude than the standard case with ν = 10 −6 v A R c . As shown in figure 10, the δṪ e component of each harmonics is compared for the linear and nonlinear phases for ν = 10 −6 v A R c and ν = 10 −8 v A R c , respectively. The poloidal structures have been modified by vortical structure with a weak dissipation value, as shown in figures 10(c) and (d), and the amplitude of δṪ e of n = 2 harmonics is significant, as shown  in figures 10(ii). Therefore, the higher viscosity coefficient is preferable to be consistent with experiments. This difference between the experiment and our simulation might be attributed to micro-turbulence. The effect of micro-turbulence on the Alfvén eigenmode evolution has been theoretically predicted to be larger than that of collisions for present tokamaks and ITER [29][30][31]. In the experiments, micro-turbulence may work as an effective collision or dissipation, and may affect the time evolution of OFM and nonlinearly generated harmonics [32].

Discussion and summary
We performed kinetic-MHD hybrid simulations to investigate the nonlinear evolution of OFM with waveform distortion in tokamak plasmas. The waveform distortion is reproduced for the first time in the nonlinear phase when the amplitude of n = 2 harmonics approaches close to that of n = 1 harmonics. It was found that the linear growth rate and mode frequency of the n = 2 harmonics are twice those of the n = 1 OFM, where the MHD nonlinearity is responsible for the generation of the n = 2 harmonics. The δṪ e signals clearly show that the emergence of waveform distortion is caused primarily by the superposition of the n = 2 harmonics on the fundamental n = 1 harmonics of OFM. We emphasize that this is a new nonlinear mechanism of waveform distortion, which differs from previous work [20]. When we looked at the time evolution of the radial profile of δṪ e , the two types of waveform with the 'falling distortion' and 'rising distortion' were found at the different regions. In particular, the 'rising distortion' is identical to the experimental observations, while the 'falling distortion' occurs closer to the plasma core and is not found in experiments. In order to understand the difference between the two types of waveform distortion, the phase difference between the n = 1 and n = 2 harmonics was analyzed using Lissajous figures, where the component δṪ e /δṪ env e were used to extract the phase as shown in figure 6. The specific Lissajous curves indicate wave couplings between the n = 1 and n = 2 harmonics with phase-lock ∼π and ∼0 for 'rising distortion' and 'falling distortion', respectively.
Through the analysis of the nonlinear components, we found that −v n=2 · ∇T eq term plays an important role in δṪ e fluctuations. For tokamak, the radial component is dominant rather than θ and ϕ components. The strong shearing structure of v rad,n=2 components at r/a ∼ 0.4 can explain the two waveform distortion types. Moreover, it should be noted that the larger viscosity is needed to reproduce the waveform distortion in the simulation than that in the experiments. This suggests that the effective viscosity in the experiments is enhanced from the collisional viscosity by some physical mechanism such as micro-turbulence.
As we mentioned in the introduction, the RWM and ELM are the biggest concern in achieving high-beta plasma. For ELMs triggered by OFMs [14,15], it is interpreted that the redistribution of energetic particles due to OFM is the main reason for the excitation of ELMs. Nonetheless, in our simulations, energetic particle redistribution does not reach the outside of plasma (r/a > 0.8), which may be related to the low mode amplitude of OFM in current cases. It would be interesting to investigate the relationship between OFM and other MHD modes by simulating the excitation of RWM or ELM following the OFM.