Multiple interactions between fishbone instabilities and internal transport barriers in EAST plasmas

Fishbone instabilities and internal transport barriers (ITBs) are frequently and sequentially observed in tokamak plasmas. Recently, the relationship between fishbone instabilities and ITBs was numerically studied, mainly on the basis of experimental results (Liu et al 2020 Nucl. Fusion 60 122001). It was identified that a radial electric field can be generated by the fishbone itself, which may act as a trigger for ITB formation. To gain a deeper understanding of this subject, in this work we further demonstrate the multiple interactions between fishbone instability and ITBs in Experimental Advanced Superconducting Tokamak (EAST) experiments (discharge #56933) using the hybrid kinetic-MHD code M3D-K. In multiple-n simulations, it is found that a zonal electric field can be induced in the nonlinear fishbone stage, leading to a relatively large E × B zonal flow that is sufficient to suppress the dominant microinstability before ITB formation; this should account for ITB triggering. After the ITB is triggered, the equilibrium pressure gradient increases and fast ions from the neutral beam injection accumulate in the ITB region. Linear simulations are performed to analyze the effect of ITB formation on fishbone instability. It is shown that due to the change of the pressure gradient during ITB expansion, the change in the bootstrap current density profile modifies the q-profile and then stabilizes the fishbone mode. Additionally, the accumulation of the fast ions leads to a broadening of fast ion distribution around the ITB region, which also has a stabilizing effect on the fishbone mode.


Introduction
Achieving steady-state operation is crucial for high performance in advanced tokamak scenarios, such as ITER [1], JT60-SA [2] and CFETR [3]. Experiments have demonstrated that the spontaneous formation of an internal transport barrier (ITB) can be used as a candidate in advanced scenarios due to the suppression of turbulent transport in the plasma core [4,5]. Such a region of enhanced confinement is associated with a high bootstrap current fraction [6][7][8], core-flat or shear-reversed safety factor q-profiles [9][10][11], strong E × B shear flow [12,13] and suppression of turbulence [14,15]. Generally, the formation of an ITB is often accompanied by magnetohydrodynamic (MHD) instability, acting as a trigger for ITB formation, such as m/n = 1/1 fishbone activities [16] and double tearing modes [17,18]. In particular, a strong relationship between fishbone instability driven by energetic particles and ITB formation has been observed in both EAST [19,20] and ASDEX Upgrade [16]. In these discharges, onset of the ITB occurs immediately after the fishbone oscillation and is radially located around q min , accompanied by a decrease in turbulent transport levels [16]. A probable explanation has been proposed, namely that the central flat q-profile can be clamped by the fishbone instability [20,21], which is considered to be a necessary prerequisite for ITB formation in experiments [9,10,22]. Such a clamping effect of the current profile [16,[23][24][25], attributed to the magnetic reconnection driven by the fishbone, can redistribute the central current as long as q drops below 1 and then clamps the q-profile with q 0 ≈ 1, which favors ITB formation. Although an interdependent relationship between fishbone instabilities and ITBs is universally accepted, further studies are required to identify the mechanism behind this interaction.
Typically, the start of ITB formation is characterized by an increase in the local density gradient or temperature gradient above the previous level. During this period, the off-axis bootstrap current related to the plasma pressure gradient is increasing, and can significantly modify the current density profile that determines the safety factor profile [6,16]. In addition, according to experimental results, the beam deposition profile is broadened when an ITB is sustained for a period of time [26]. A fishbone instability is a MHD instability driven by energetic particles [27][28][29]; from the viewpoint of macroscale instabilities, any change in MHD equilibria or kinetic parameters during ITB evolution is bound to influence fishbone linear instability, such as a change in the pressure gradient or energetic particle distribution fraction. From the perspective of microinstabilities, the anomalous transport level is a determining triggering factor of ITB formation. Anomalous transport can be suppressed by zonal flow generated by turbulence itself, which has been studied extensively [12,30,31]. Further, it has been confirmed that a oscillating zonal flow can be induced by mutual interaction between energetic particles and MHD modes [32], such as fishbone instabilities [27] and Alfvén modes [33]. Previous experimental works have proposed that a radial electric field is induced by the fishbone and appears near the foot point of the ITB, then a corresponding E × B shear flow is generated [15,16]. Once the E × B shear flow is strong enough to suppress the dominant microinstability, the transport levels decrease [13,14] and an ITB is formed. However, this possibility has not yet been systematically validated. Motivated by the analyses above, several possible interactions between fishbone instabilities and ITBs in EAST are numerically studied in this work (figure 1).
In this paper, it is found that a fishbone becomes stable with increasing bootstrap current caused by equilibrium variation during ITB expansion, which provides an explanation for the disappearance of the fishbone in the later ITB phase. It is numerically verified that the broadening profile of the energetic ion distribution also contributes the stability of the fishbone. To analyze the radial electric field of the n = 0 component driven by mode coupling during fishbone activities, a multiple-n simulation is performed. It is shown that E × B zonal flow can be induced during nonlinear fishbone evolution, and the maximum E × B shear flow appears near the q = 1 rational surface where the ITB is to be formed. Most importantly, the E × B shearing rate is always greater than the growth rate of the dominant microinstability during nonlinear evolution of the fishbone, which is considered to be an important trigger of ITB formation. The remaining parts of this paper are organized as follows. Section 2 describes the experimental observations in EAST, the model of the M3D-K code and the equilibria used in this work. In section 3, linear simulations show the effect of ITB formation on fishbone instability, including the impact of changes in the q-profile and energetic particle distribution. To study the triggering effect of the radial electric field and the shear flows generated by the fishbone on ITB formation, nonlinear simulations are carried out in section 4. Section 5 contains conclusions and discussion concerning this work.

Experimental observation and simulation setup
A typical high-β N discharge (EAST #56933) with a strong interaction between fishbone instability and the ITB is shown in figure 2. It was illustrated in detail in experimental work for this shot [34]. The start time of the ITB was chosen at EAST discharge #56933. Time evolution of (a) NBI heating power and 4.6 GHz LHW heating power, (b) normalized beta (β N ), (c) multichannel ECE signals and (d) the wavelet spectra (analysis with soft x-rays). (e) The poloidal structure of the fishbone at t = 3.65 s from SX tomography. t = 3.59 s due to the increased core gradient of T e within ρ < 0.5, as observed in the multichannel electron cyclotron emission (ECE) signals of figure 2(c), where ρ is the square root of the normalized toroidal magnetic flux. Note that the ECE power displayed here is only used to analyze the relative change of T e rather than the absolute value of T e , as explained in previous work [34]. As shown in figure 3 of Gao et al [35], with auxiliary lower-hybrid wave (LHW) and neutral-beam injection (NBI) heating, it was specified that strong and broad ITBs were formed in n e , T e and T i profiles at t = 3.59 s. Then, Figure 3. Growth rate and frequency at ρ = 0.36 before (t = 3.55 s) and after (t = 3.75 s) ITB formation for EAST #56933, calculated by the transport code TGLF. Here k y is the normalized poloidal wave number.
the ITB was fully developed up to t = 3.75 s, located near ρ = 0.25-0.5. Notably, after the ITB was fully developed, this typical discharge #56933 achieved the largest value of the triple product n i0 T i0 τ E ∼ 1.0 × 10 19 m −3 keV s for EAST so far. After t = 3.75 s, the ITB was gradually degraded. Before ITB formation, fishbone MHD activity was observed clearly in the spectra with a frequency of 3-5 kHz at t = 3.58 s, which is near the start time of ITB formation, as shown in figure 2(d). Then, fishbone activity disappeared before ITB decay at t = 3.68 s. In figure 2(e), the poloidal structure of the fishbone at t = 3.65 s, given by SX tomography, shows an obvious twisting mode structure pattern, and its poloidal mode number is clearly identified as m = 1. To analyze the microinstability in shot #56933, we use the transport code TGLF [36,37] to depict the k spectrum of the linear growth rate and frequency of the most unstable modes at the core gradient area ρ = 0.36 before (t = 3.55 s) and after (t = 3.75 s) ITB formation, as shown in figure 3. The normalized growth rate, real frequency and poloidal wave number are γ = γ(a/c s ), ω = ω(a/c s ) and k y = k θ ρ s , respectively, where c s = T e /m i , ρ s = c s /Ω s , Ω s = eB/m i c and a is the circular equivalent minor radius of the last closed flux surface. In figure 3(b), the positive (negative) value of frequency (ω ) indicates that the mode rotates along the electron (ion) diamagnetic drift direction. Figure 3 shows that the dominant unstable modes in both cases are low-k (k y < 1) modes, rotating along the ion diamagnetic drift direction, which should be the ion temperature gradient (ITG) modes [38,39] before and after ITB formation. In particular, the growth rate of the ITG mode decreases after ITB formation in figure 3(a), contributing to the improved confinement and ITB formation. However, although the reduction in the ITG growth rate was considered to be the leading factor in ITB formation in this discharge [35], the reasons for this reduction have not been further discussed. This point will be analyzed below. It has also been suspected that the fishbone's disappearance is caused by the change in core q-profile in this discharge [34,40]. However, measurement of the q-profile has a lot of uncertainty in experiments, especially in the absence of MHD activities as the reference to modify the q-profile after ITB formation. In this work, therefore, a simplified model is used to analyze how the q-profile changes during ITB formation and the impact of such change on fishbone stability. Moreover, although fishbone instability has generally been accepted as a trigger for ITB formation in experiments, the mechanism behind this still requires detailed discussion. Our preliminary working hypothesis is that the increase in bootstrap current, generated by the steepening pressure gradient in the ITB region, changes the plasma total current and therefore the q-profile, thus altering fishbone stability. In turn, the zonal flow induced in the nonlinear fishbone stage may contribute to the formation of an ITB. In this work, we apply the global kinetic-MHD hybrid nonlinear code M3D-K [27,41] to verify our hypotheses based upon the experimental equilibria in shot #56933. For the MHD part in code M3D-K, the dissipative MHD model consists of the full resistive MHD equations and the additional dissipative terms including viscosity and heat conduction. The M3D-K code has been widely used to analyze the linear and nonlinear behaviors of energetic particle modes and Alfvén eigenmodes in tokamak plasmas [42][43][44]. Figure 4 shows the equilibrium profiles from EAST discharge #56933, including plasma pressure, density n e , electronic temperature T e , ion temperature T i and safety factor q-profile. The blue curves denote the equilibrium profiles before ITB formation at t = 3.55 s, approximately 0.03 s prior to appearance of the fishbone. The red curves represent the equilibria after the disappearance of the fishbone in ITB formation later at stage at t = 3.75 s. For this discharge, the profiles of T e , n e and T i were measured by Thomson scattering (TS) and x-ray crystal spectrometer (XCS) diagnostics and reflectometry. It is important to note that the highest time resolution of the TS diagnostics is 100 ms. More details of diagnostic methods for these profiles can be found in Yang et al's experimental work [34]. In figure 4(e), the safety factor q-profile at t = 3.55 s is obtained from the EFIT [45] reconstruction combined with the Faraday rotation angle measured by the POINT (polarimeter-interferometer) [46] diagnostics. Several experimental MHD observations, such as the previous sawtooth oscillation and the area of the fishbone, are also considered in q-profile corrections. Note that the minimum value of the weak reversed shear safety factor q is slightly below 1. It can be observed that the q-profile is flat in the plasma core with q 0 ≈ 1; such a core-flat q-profile with integer q rational surface has been confirmed to be favorable for ITB triggering and maintenance by a series of experimental reports, such as on ASDEX Upgrade [47], JET [9,10,22], JT-60U [48] and DIII-D [49] devices, and even on a stellarator [50]. The main parameters are as follows: inverse aspect ratio ε = a/R = 0.237, triangularity δ = 0.43, elongation κ = 1.6, on-axis magnetic field B 0 = 1.58 T, electron central density . Spatial equilibrium profiles of (a) normalized plasma pressure, (b) normalized density, (c) normalized electron temperature, (d) normalized ion temperature and (e) safety factor q-profile. The blue and red curves denote the equilibrium profiles before (t = 3.55 s) and after (t = 3.75 s) ITB formation, respectively. In (e), the blue curve and blue dotted line represent the safety factor q = 1 and its corresponding radial location, respectively, and the shaded area indicates the ITB region. Note that the q-profile at t = 3.55 s is reconstructed from the equilibrium fitting code EFIT.
Alfvén frequency ω A = εv A /R 0 , central total beta value β tot,0 = 3.3% and the ratio of central energetic ion beta to the central total plasma beta is β h,0 /β tot,0 = 0.2. The plasma viscosity is chosen to be 10 −6 . A static equilibrium is calculated by the equilibrium code VMEC [51] based on the above profiles and parameters, and is used as the initial condition for the M3D-K code. It is worth pointing out that there is a degree of freedom in the choice of inputting equilibrium profiles to find the plasma equilibrium in the VMEC code. Two commonly used approaches are q-solver and J-solver [52][53][54], which obtain the plasma equilibrium set by specifying q, P or J, P profiles, respectively. Here P, q and J are plasma pressure and safety factor and total toroidal current density, respectively.
For the NBI heating system in discharge #56933 there are two beam injectors, one for co-current injection and another for counter-current injection. Each injector consists of two ion sources. The NBI power was applied at 2.4 s and then increased stepwise by presetting the injection time sequence of the four ion sources. At t = 3.55 s, the NBI power was 1.1 MW [34] and the corresponding energy of the co-current beam injection was 60 keV [55]. Based on these details, the properties of the energetic ions adopted in our simulations are given below. The isotropic distribution of energetic ions is slowing down in velocity space and peaks in pitch angle, and its form is given as follows: where C is a normalization factor, H is the step function, v 0 is the neutral beam injection speed, v c is the critical velocity, Λ ≡ μB 0 /E is the pitch angle and μ is the magnetic moment. The central pitch angle value of energetic ions is selected as Λ 0 = 1.0 and the width of the pitch angle is ΔΛ = 0.5. Notably, the energetic ion distribution includes a significant fraction of passing particles when Λ 0 = 0 and a major fraction of trapped particles when Λ 0 = 1. Ψ is the normalized poloidal magnetic flux Ψ averaged over the energetic ion orbit and ΔΨ = 0.3. The normalized gyroradius and speed of energetic ions are chosen as ρ h /a = 0.065 and v 0 /v A = 0.686, respectively.

The stabilizing effect of ITB expansion on the fishbone
In this section, based upon the experimental equilibrium profiles before ITB formation at t = 3.55 s, the numerical results for the m/n = 1/1 fishbone are looked at first. Figure 5 shows that the linear mode structure (U) is located inside the q = 1 rational surface and has a twisted character, agreeing well with the experimental measurements in figure 2(e). Here U is the velocity stream function and it is associated with the incompressible part of the plasma velocity by the equation Further, even after minor adjustment of the q-profile in this work, the frequency behavior of the fishbone is similar to that in a previous simulation study by Liu et al [19], which was confirmed to be consistent with experimental observations [34]. The successful reproduction of a fishbone in this work forms the basis for the further analyses of fishbone instability during ITB expansion.
It has been reported that the large bootstrap current generated by the steep plasma pressure gradient in the ITB period allows steady state operation [16]. During ITB expansion, due to the change in the bootstrap current density profile caused by variation of the plasma pressure gradient, the total current density profile can be altered. Then, the change in q-profile determined by the total current density will affect fishbone stability [7,56]. In this section we now analyze the effects of q-profile variation due to the change in bootstrap current on fishbone instability during ITB formation.
Before ITB formation stage at t = 3.55 s, the occurrence of the m/n = 1/1 fishbone clamps the central flat q-profile locally, keeping the central q-value q(0) in the vicinity of 1 [21,34], as shown in figure 4(e). For the later ITB formation stage at t = 3.75 s, a reduced model is developed as follows to obtain the q-profile due to uncertainty in the q measurement in experiments. First, according to the experimental q-profile at t = 3.55 s, the total toroidal current density profile at t = 3.55 s can be obtained with q-solver, as shown by the solid curve in figure 6. Then, we make a simplifying assumption that the total toroidal current consists of two parts: the bootstrap current and the driven current. Here, the driven current in our simulation is mainly composed of three parts, ohmic current, NBI driven current and LHW driven current, which is roughly consistent with the current compositions in EAST #56933 [34]. In addition, by analyzing the current profiles during the fishbone in EAST #71320 [21], it is found that the ohmic current accounts for 80% of the total driven current, meaning that the variation in the driven current is mainly dominated by the ohmic current. Notably, discharge #71320 applied similar heating powers to #56933 to confirm the current profiles in #56933 [34]. In this work, the overall plasma toroidal current I tot is kept constant. In addition, according to the estimation of the neoclassical/classical current diffusion time in EAST #56933 in previous work [21], the current diffusion time of the neoclassical bootstrap current is much shorter than that of the classical driven current. That is, compared with the driven current, the bootstrap current has a shorter response time for the change in equilibrium profile. Hence, in our reduced model, the shape of the driven current density profile J d is considered to be almost unchanged during a sustained ITB. Note that the driven current I d decreases with increasing bootstrap current I bs due to the constant toroidal current I tot . The decreased driven current I d has been confirmed by the decreased ohmic heating power in #56933. Therefore, evolution of the total toroidal current profile is mainly dependent on the bootstrap current, which is related to the plasma equilibrium profiles and will be given below. The profile of the driven current density is obtained by subtracting the bootstrap current density from the total toroidal current density. Finally, using the obtained total toroidal current density profile at t = 3.75 s, the q-profile at this moment can be predicted by J-solver. In the following, the detailed calculation procedure for this model is presented.
According to previous work [8], the bootstrap current density profile is given by where J bs is the bootstrap current density, ε is the inverse aspect ratio, p e is the electronic pressure, B θ is the poloidal component of the magnetic field, ψ is the poloidal magnetic flux and n and T are the plasma density and temperature, respectively.
The distribution of bootstrap current density at t = 3.55 s is calculated from equation (3), as shown by the dashed curve in figure 6. It is found that the bootstrap current density profile peaks at around ρ = 0.3. According to previous experiments in EAST discharges [20,21,57], the bootstrap current fraction at t = 3.55 s is set to be I bs /I tot = J bs r dr/ J tot r dr = 27% in our study. Therefore, the distribution of the driven current density is shown by the dotted curve in figure 6.
Similarly, the bootstrap current density at 3.75 s can also be obtained via equation (3), and its profile is represented by the red curve in figure 7(a). As expected from the results of figure 7(a), compared with the time at 3.55 s, the steeper pressure gradient leads to an increase in the bootstrap current density, peaking at t = 3.75 s. As mentioned above, the fraction of bootstrap current is premised to be 27% at t = 3.55 s and the value of the total plasma toroidal current remains unchanged during ITB formation. Thus, the fraction of bootstrap current at t = 3.75 s is calculated to be 45% based on the relative change of bootstrap current density at t = 3.55 s and at t = 3.75 s. Then, combining the driven current density profile J d , whose amplitude decreases, the total toroidal current density profile at t = 3.75 s is obtained, as shown in red curve of figure 7(b). It is shown that with the equilibrium calculation based on a reduced model during ITB expansion, the total toroidal current density decreases in the plasma core and increases in the ITB position near ρ = 0.4. Then the q-profile at t = 3.75 s is obtained by J-solver, as indicated by the red curve in figure 7(c). From t = 3.55 s to t = 3.75 s, it is found that the central safety factor q(0) is elevated while the value of q near ρ = 0.6 slightly decreases. Notably, at t = 3.75 s, the q = 1 rational surface disappears and q min = 1.04. As the bootstrap current fraction increases during ITB formation, the reversed magnetic shear increases. Up to now, the q-profile after the fishbone disappearance at t = 3.75 s has been calculated by the above analyses with our reduced model. Subsequently, the equilibrium profiles at t = 3.75 s were employed to analyze the mode instability through M3D-K, and the results are presented in figure 8. It is shown that the onedimensional mode structure at t = 3.75 s is located inside the q min rational surface and does not change much compared with that at t = 3.55 s in figure 8(a). Further, the linear numerical results of mode frequency do not differ significantly between these two moments in figure 8(b), which implies that the unstable mode at t = 3.75 s is still a fishbone. However, the mode growth rate reduces by half from t = 3.55 s to t = 3.75 s. Such a change in fishbone stability, caused by equilibrium variation during ITB growth, is consistent with the phenomenon in EAST #56933, giving an explanation for the disappearance of the fishbone before 3.75 s in the experiment. From the above simulation results, the physical process can be understood as follows: with the steepening plasma pressure gradient during ITB development, the bootstrap current fraction increases, leading to the increase in reversed magnetic shear and the complete disappearance of the q = 1 rational surface. Then the fishbone instability becomes stable due to the disappearance of the q = 1 rational surface that is crucial for classical fishbone instabilities, which may cause the consequent degradation of the ITB. In addition, since the integer q rational surface plays a significant role in ITB formation and maintenance, as mentioned in section 2, the disappearance of the q = 1 surface may also be a direct cause of ITB decay. In fact, experiments have proposed that the disappearance of the fishbone and the consequent decay of the ITB may be caused by the change in the q-profile, for example further change to reversed shear [34,40]. Due to the lack of accurate measurements of the q-profile after disappearance of the fishbone, no further analyses were conducted. In this study, on the basis of the q-profile before ITB formation at t = 3.55 s and calculation of the bootstrap current based on the experimental pressure profiles before and after ITB formation, the q-profile after ITB expansion at t = 3.75 s is obtained with some simplifying assumptions. Our linear simulation results show that the variation in the q-profile caused by the change of bootstrap current during ITB expansion has a stabilizing effect on the fishbone, supporting the prediction that disappearance of the fishbone may be caused by the change in the q-profile during ITB expansion. It is worth mentioning that, after the sensitivity analyses, our results are robust within the highest bootstrap current fraction (i.e. 50%) that EAST can achieve at present. Besides, considering there is a certain gap between our reduced calculation process and the real experimental evolution, the analyses are mainly qualitative rather than quantitative.
According to experiment [26], the ITB sustainment phase is not exactly a steady state. Not only does the q-profile change slowly in time but the beam deposition profile also gradually broadens during this period. Figure 9(a) presents the energetic ion density distribution calculated with the Monte Carlo code NUBEAM [58] based on the experimental data in shot #56933. During ITB formation from 3.55 s to 3.75 s, it is found that the deposition position of energetic ions is around ρ = 0.2-0.4, which is near the region of ITB formation, leading to a broader energetic ion distribution within the q = 1 rational surface. Based on this result, the effect of the parameter Δψ on fishbone instability is studied, where Δψ is the radial width of the energetic ion distribution in equation (1). The value of Δψ ranges from 0.2 to 0.8 in our simulation. Here, two examples of Δψ at 0.3 and 0.8 are shown in figure 9(b). To clearly study the effect of the broadening energetic ion distribution on the fishbone instabilities, only Δψ is changed in these simulations, with the other equilibrium parameters being the same as those at t = 3.55 s.
The linear growth rate and mode frequency versus different radial widths of the energetic ion distribution Δψ are displayed in figure 9(b). It is found that the growth rate first increases from Δψ = 0.2 to Δψ = 0.3 and then rapidly decreases with increasing Δψ. This result implies that the fishbone can be stabilized by broadening the energetic ion distribution, which may also serve as one of the reasons for the experimental disappearance of the fishbone in later ITB stages. The relevant physical mechanism behind this process can be understood as follows. As the ITB is formed, the turbulent transport level decreases and the energetic ions accumulate around the ITB region, which leads to broadening of the energetic ion distribution within the q = 1 rational surface. Since the fishbone instability is mainly driven by d f /dP φ [27] of energetic particles located within the q = 1 surface [27,59], where f denotes the energetic ion distribution and P φ is the toroidal angular momentum, the mode instability drive becomes weaker as the energetic ion distribution broadens. The fishbone then eventually becomes stable in this process. Further, the mode frequency slightly decreases as Δψ increases, which is probably due to the resonant particle changes resulting from broadening of the energetic ion distribution.

The triggering effect of nonlinear fishbone evolution on ITB formation
To date, the above simulations only includes energetic particle nonlinearity, with MHD nonlinearity being neglected. In this section, we mainly focus on the nonlinear dynamics of the fishbone, which is important for understanding the experimental observations as well as the physical mechanisms behind these observations. Consequently, both particle nonlinearity and MHD nonlinearity are included in the following simulations to study the impact of nonlinear fishbone evolution on ITB formation. We apply the same equilibrium profiles and parameters as described in section 2 at t = 3.55 s, and the discrete toroidal modes n = 0-4 are included. Here, the n = 0 component is mainly self-generated by nonlinear coupling of the n > 0 modes to themselves. Furthermore, the radial electric field of the n = 0 component is the main focus. Figure 10 contrasts the time evolutions of kinetic energy in a single-n simulation (dotted curve) and a multiple-n simulation (solid curves). For the multi-toroidal mode case, it is shown that the n = 1 mode is still the dominant mode, indicating that the n = 0 component is mainly generated by mode coupling of the n = 1 modes. It can be seen that the mode kinetic energy amplitude of n = 0 is comparable with the n = 1 component. In addition, the MHD nonlinearity does not enhance the initial mode growth rate of the n = 1 component but reduces the initial saturation level of the n = 1 mode by more than half. These results are consistent with the classical (1, 1) fishbone results [27].
It is well known that E × B zonal flow plays an important role in turbulence regulation and ITB formation as it contributes to transport reduction [12,[60][61][62][63]. Further, it has been reported that a zonal electric field can be induced by the mutual interaction between energetic particles and the MHD mode [32], such as fishbone instabilities. Therefore it is necessary to study the link between E × B zonal flow induced in the nonlinear fishbone stage and ITB formation. In the following analyses, n is the toroidal mode number in the multiple-n case if not specially noted. In addition, according to figure 10, t = 700τ A and t = 1950τ A are defined as the moments when the fishbone enters the initial saturation stage and enters the later saturation stage, respectively. Figures 11(a) and (b) show the n = 0 component of |E r | at t = 700τ A and t = 1950τ A , respectively, where |E r | = |e r · ∇Φ| is the radial electric field amplitude and Φ is the electric potential. It is shown that the zonal electric field is localized in the plasma core in the initial saturation phase in figure 11(a). In the later saturation phase in figure 11(b), the electric field moves outwards and appears near the q = 1 rational surface, located inside the ITB formation region with ρ = 0.25-0.5. This numerical result implies that the zonal electric field induced in the later saturation phase of the fishbone can appear near the q = 1 surface, where the ITB is to be formed. According to equation (2), the electric field structure has a close relationship with the mode structure. To better understand the mechanisms behind the time evolution of the zonal electric field, figures 11(c) and (d) show the mode structures of the n = 1 component in the initial saturation phase and in the later saturation phase, respectively. In figure 11(c), the mode structure of the fishbone in multiple-n simulation in the initial saturation phase shows a similar result to that in single-n simulation in the linear phase (figure 5), located inside the q = 1 surface. Since the n = 0 component is mainly generated by mode coupling of the n = 1 mode to itself, the zonal structure is related to the n = 1 mode structure. Thus, the zonal electric field induced in the initial saturation phase of the fishbone is localized in the plasma core in figure 11(a). However, the mode structure moves radially outward and shrinks around the q = 1 surface in the later saturation phase, as shown in figure 11(d). Such an evolution of mode structure during the nonlinear fishbone process can be explained by energetic ion redistribution induced by nonlinear wave-particle trapping effects [27,64]. Figure 11(e) shows the energetic ion radial distribution at t = 0, 700 and 1950τ A . In the initial saturation phase at t = 700τ A , it is found that the region of energetic ion redistribution is mainly located inside the q = 1 surface. Compared with the redistribution at t = 700τ A , the level of energetic ion redistribution is stronger at t = 1950τ A , and the flattening region of the energetic ion distribution extends beyond the q = 1 surface. This result indicates that there are considerable energetic ions to be redistributed outside the q = 1 surface in the later saturation phase, leading to change of the wave-particle resonance region and therefore the fishbone mode structure. Then, the zonal electric field near the q = 1 surface, located within the ITB formation region (ρ = 0.25-0.5), can be induced in the later saturation phase of the fishbone, as shown in figure 11(b), and thus E × B zonal flow will also appear nearby. Next, we will carry out a quantitative analysis of E × B zonal flow.
Although the above results suggest that a radial zonal electric field can be driven by fishbone activity near the q = 1 surface, the E × B shear flow caused by this electric field is the accepted mechanism for suppressing the microinstability [63,65]. The E × B shearing rate for turbulence suppression is defined as [63] where R and r are the major and minor radii, respectively, B θ and B φ are the poloidal and toroidal magnetic components of the magnetic field, respectively, and E r is the radial electric field that we obtained above. According to the associated report [66], when the E × B shearing rate significantly exceeds γ max in the core region, suppression of microinstability may occur, leading to a reduced transport level and improved confinement, and then an ITB will be formed. Here, γ max is the maximum linear growth rate of the dominant microinstability in the plasma. Even when γ E×B ≈ γ max , in particular, nonlinear gyro-Landau fluid simulations have also proved complete stability of turbulence [67,68]. This simplified criterion has thus been widely used for quantitative comparison between experiment and theory [69]. In order to quantify the comparison, the shearing rate γ E×B is calculated according to equation (5). It is found that during nonlinear fishbone evolution the maximum γ E×B is located around the q = 1 surface with ρ = 0.4, near the region where the ITB is to be formed and close to the position with ρ = 0.36 where the linear growth rate of microinstability is calculated by TGLF in figure 3. To illustrate this, two examples of the γ E×B structure at t = 700τ A and at t = 1950τ A are given in figures 12(a) and (b), respectively. Then we calculate the time evolution of the γ E×B value near the q = 1 surface in figure 13. Note that the fishbone at t = 600τ A is still in the linear stage rather than the nonlinear stage. In figure 13, the shearing rate γ E×B increases significantly when the fishbone reaches the initial saturation stage at t = 700τ A compared with the linear stage at t = 600τ A , after that, γ E×B gradually decreases. Such a trend of variation of the shearing rate should be associated with the saturation amplitude evolution of the fishbone mode in figure 10. The red dotted line in figure 13 represents the maximum linear growth rate γ ITG, max of the ITG instability at ρ = 0.36 and t = 3.55 s on EAST #56933, obtained by TGLF in figure 3. It can be seen that γ E×B is consistently greater than γ ITG, max during nonlinear fishbone evolution, and γ E×B can even be several times larger than γ ITG, max in the initial saturation stage of the fishbone. According to the criterion of turbulence suppression by E × B shear flow that γ E×B > γ max ,  the E × B zonal flow induced in the nonlinear fishbone stage is sufficient to suppress the ITG instability in shot #56933, which is required for the formation of an ITB. As long as an ITB is formed, the microinstability is further suppressed, leading to the decreased linear growth rate of the ITG mode after ITB formation at t = 3.75 s in figure 3. Actually, fishbone activity has been suggested to act as a trigger for ITB formation in ASDEX Upgrade in a previous study [16]. One possible mechanism suggested by the authors was that the redistribution of energetic particles during fishbone oscillation causes a current across the rational surface, resulting in a sheared radial electric field with a shearing rate comparable to the growth rate of ITG modes, which may be responsible for the suppression of turbulent transport. Our numerical results visually show a support for the likelihood that the radial electric field is strongly associated with the redistribution of energetic particles during fishbone activities. Most importantly, it is quantitatively confirmed that the zonal flow induced in the nonlinear fishbone stage is indeed enough to suppress the turbulent instabilities, which is a key trigger for ITB formation in EAST #56933.

Conclusion and discussion
In this work, linear and nonlinear simulations of interdependent relationships between fishbone instability and ITBs on EAST are systematically studied with the global hybridkinetic code M3D-K. To understand the impact of ITB expansion on linear fishbone instability, two possible scenarios during ITB formation are discussed. One scenario is that the change in the bootstrap current density profile, caused by the variation of equilibrium pressure during ITB growth, will result in disappearance of the q = 1 rational surface, while the q = 1 surface is essential for both fishbone existence and ITB formation. The fishbone becomes stable during this process. The other scenario is that the reduction of the transport levels after ITB formation leads to the accumulation of energetic ions in the ITB region, which broadens the distribution of energetic ions. This process is also supported by the experimental data. The numerical results show that the fishbone gradually becomes stable with broadening energetic ion distribution. Both the above scenarios can account for vanishing of the fishbone during ITB expansion in experiments. In the nonlinear multiple-n simulation, the MHD nonlinear effects reduce the initial saturation level of the fishbone. A localized electric field of the n = 0 component induced in the nonlinear fishbone stage appears near the region where the ITB is to be formed, indicating that E × B zonal flow can also be induced around here. Moreover, the E × B shearing rate γ E×B is always higher than γ ITG, max for shot #56933 in our simulations, playing a dominant role in the triggering and growth of the ITB, where γ ITG, max is the maximum linear growth rate of the dominant microinstabilities. As soon as the ITB is formed, the increased pressure gradient provides a positive feedback in the increase of γ E×B to further suppress the microinstabilities and radial transport during ITB growth, giving a reasonable explanation for the decreased γ ITG, max after full development of the ITB in the experiment.
Therefore, the results of the multiple interactions between fishbone instabilities and ITBs in EAST plasmas shown in our study can be summarized as follows. After fishbone activity appears, E × B zonal flow can be induced during nonlinear fishbone evolution and eventually appears near the region where the ITB is to be formed. In this process, the shearing rate of E × B zonal flow is consistently greater than the upper limit of the microinstability drive, which is the most likely candidate for confinement transition and acts as a key trigger for the ITB. Once the ITB is formed, the pressure gradient increases in the ITB region, providing a positive feedback to further increase E × B zonal flow, leading to decreased transport levels in this region, which contributes to further development of the ITB. Besides, during ITB expansion, both the increased reversed shear of the q-profile and the broadening energetic ion distribution will stabilize the fishbone.
Our simulations not only provide explanations for the experimental phenomenon of fishbone disappearance during ITB expansion, but also further prove that E × B zonal flow can be induced during nonlinear evolution of fishbone instability. In particular, our simulation results also quantitatively confirm that the E × B zonal flow induced in the nonlinear fishbone stage is sufficient to suppress the dominant microinstability in EAST #56933; we thus interpret this result as ITB formation. But in fact, even without considering the influence of the mode saturation level, the obtained numerical results for the E × B shearing rate are based on slight adjustment of the experimental equilibrium profiles. Additionally, when carrying out the quantitative comparison, certain errors inevitably exist between the E × B shearing rate calculated by M3D code and the linear growth rate of the microinstability calculated by the TGLF code. Even then, using typical experimental equilibrium profiles and parameters, the E × B shearing rate obtained in our simulations is comparable to or even several times larger than the maximum linear growth rate of the dominant microinstability in EAST #56933. Compared with our previous work [19], there are several obvious extensions in this work, such as analysis of the bootstrap current during ITB expansion and of the E × B zonal flow induced in the nonlinear fishbone stage. It should be noted that the particle transport process and the time evolution of density perturbations are not considered in this work; these can affect E × B zonal flow and will be the focus of future self-consistent nonlinear simulations.