# A self-sustaining mechanism for Internal Transport Barrier formation in HL-2A tokamak plasmas

W H Lin, J Garcia, J Q Li, S Mazzi, Z J Li, X X He, X Yu

► To cite this version:

W H Lin, J Garcia, J Q Li, S Mazzi, Z J Li, et al.. A self-sustaining mechanism for Internal Transport Barrier formation in HL-2A tokamak plasmas. Nuclear Fusion, 2023, 63, pp.126048. 10.1088/1741-4326/ad0508. cea-04823534

HAL Id: cea-04823534

https://cea.hal.science/cea-04823534v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A self-sustaining mechanism for Internal Transport Barrier formation in HL-2A tokamak plasmas

W. H. Lin¹, J. Garcia², J. Q. Li¹,ᵇ, S. Mazzi², Z. J. Li¹, X. X. He¹, X. Yu¹

¹Southwestern Institute of Physics, Chengdu 610041, China.

²CEA, IRFM, F-13108 Saint Paul-lez-Durance, France.

ᵇAuthor to whom correspondence should be addressed: lijq@swip.ac.cn

(Dated: November 27, 2023)

The formation of Internal Transport Barrier (ITB) is studied in HL-2A plasmas by means of nonlinear gyrokinetic simulations. A new paradigm for the ITB formation is proposed in which different physics mechanisms play a different role depending on the ITB formation stage. In the early stage, fast ions, introduced by Neutral Beam Injection (NBI) ion system, are found to stabilize the thermal-ion-driven instability by dilution, thus reducing the ion heat fluxes and finally triggering the ITB. Such dilution effects, however, play a minor role after the ITB is triggered as electromagnetic effects are dominant in the presence of established high pressure gradients. We define the concept of ITB self-sustainment, as the low turbulence levels found within the fully formed ITB are consequences of large scale zonal flows, which in turn are fed by a non-linear interplay with large scale high frequency electromagnetic perturbations destabilized by the ITB itself.

## I. INTRODUCTION

The final goal of magnetic confinement devices is to confine plasmas of high temperature and density for sufficiently long time in order to produce economically advantageous fusion energy. Confined plasmas can be severely degraded by the outward energy transport driven by micro-instabilities such as the Ion-Temperature-Gradient (ITG) mode [1]. Therefore, a credible path towards reliable energy fusion production must rely on mechanisms controlling such an energy transport.

Plasmas with Internal transport barriers (ITB) [2], characterized by a suppression of heat transport driven by microturbulence leading to high core temperatures and densities, have been shown to provide a way to improve plasmas energy confinement in various tokamaks [3–8]. The formation and characteristics of ITB have been extensively studied. Several physical mechanisms have been put forward to explain energy transport reduction or suppression within an ITB. One of the initial mechanisms proposed was the $E \times B$ flow shear turbulence stabilization (see [2] for example), which manifests itself by breaking up turbulent eddies and reducing the amplitude and cross phase of turbulent fluctuations. In this context, negative or low magnetic shear is also known to have a synergistic effect with $E \times B$ shear on ITB formation, as it weakens the drive of some unfavourable instabilities [9] on one hand and prevents the detrimental effects brought by $E \times B$ shear [10] on the other. Other mechanisms related to the presence of highly energetic fast ions have been proposed as well. A large fraction of fast ions produced from neutral beam injection (NBI) are found to be crucial in ITB formation by their dilution effects [11], while a small minority of them could also be decisive through mechanisms such as linear resonant interaction with ITG [12] or the enhancement of α-stabilization [4].

Despite the amount of studies devoted to clarify the

physical mechanism behind the ITB formation, there are still aspects that remain unclear, e.g., whether a single physical mechanism or multiple ones are responsible for the ITB triggering and whether such mechanisms play significant roles on the ITB sustainment once it is fully formed. Clarifying these aspects is essential in order to properly evaluate whether plasmas with ITBs will be possible in future fusion reactors, for which some mechanisms, such as the $E \times B$ shearing produced by external injected torque, are known that will be less efficient.

In this work, it is shown that the triggering and sustainment of ITB rely on two different physical mechanisms depending on the ITB formation stage. Whereas the ITB triggering is found to be a consequence of the NBI fast-ion dilution, it is proposed the concept of self-sustainment of the ITB as it is the ITB itself producing the physical mechanisms that provides its sustainment. The increase of electromagnetic (EM) effects in the presence of strong ITB-generated pressure gradients reduces turbulence and transport through the onset of large scale zonal flows [14] (with toroidal number $n = 0$ and frequency $\omega = 0$), which tap energy non-linearly from large scale MagnetoHydroDynamics (MHD) fluctuations that are destabilized by the ITB itself. Meanwhile, the $E \times B$ shearing generated by the plasma rotation is not found to play a major role on the ITB formation. Such findings may pave the way for the formation of ITBs in future tokamaks as long as EM effects are dominant.

Turbulence and transport analyses are performed with state-of-the-art gyrokinetic simulations for the ITB discharge #22453 in the HL-2A tokamak [8]. In such a discharge as shown in Fig. 1, the ITB is triggered at about $t=510$ms, and stably sustained for a time window of about 250 ms. During the ITB formation, the core ion temperature has increased from 1.0 to 2.3 keV, forming a region of large $R/L_{T_i}$ ($\approx 20$) with the ITB foot located at $\rho_{tor} \approx 0.4$. Here, R is the major radius, $L_{T_i}$ the inverse logarithmic gradient of ion temperature and
<figure><img src="image_2.png" /><figcaption>FIG. 1: (a) Profiles of ion temperature $T_i$ measured by CXRS [13] at different time points in HL-2A discharge #22453. (b) The frequency spectrum of mirnov coil signal and the normalized logarithmic ion temperature gradient $R/L_{T_i}$ at $\rho_{tor} = 0.25$ in the same discharge.</figcaption></figure>

$\rho_{tor}$ the normalized square root of toroidal magnetic flux. The profiles at 510 and 650 ms, when the ITB begins to trigger and has been fully developed, respectively, are of particular interest to our study and provide the parameters set for the simulations discussed below. Shortly after the ITB triggering, the Mirnov coils detect two perturbations, a weak one at 70 kHz and a stronger one at 20 kHz in the laboratory frame, the latter being identified as a long-lived mode (LLM) in previous work [15]. LLMs, as well as fishbone (FB) instabilities [16], are both MHD modes frequently observed in HL-2A after ITB triggering. Although FBs are proposed as the key factors of the ITB formation in some tokamaks [17–19], they could hardly be related to the ITB triggering in HL-2A [20] where FBs are less observed preceding the ITB. As for LLMs, a limited amount of works [21] exist regarding their effects on ITB. In spite of the fact that n=1 LLMs often emerge near $q = 1$ rational surface in HL-2A ITB plasmas, LLMs as the non-resonant kink modes [22] do not rely on the existence of $q = 1$ surface to be destabilized. Therefore, the emphasis of this work, unlike some previous works [23–26], is placed on the dynamic interplay between ITB and LLM rather than on the rational

surfaces.

The structure of this paper is arranged as follows: after the simulation setup is addressed in section 2, the dominant instabilities in various simulation conditions are analyzed in section 3, and the stabilizing factor that is of vital importance on ITB formation is investigated by the analysis of ion heat flux in section 4. It will be shown that the full ITB formation benefited from not only the linear stabilization of dominant instabilities but also the nonlinear EM effect, which is attributed in section 5 to the onset of zonal flow through the saturation of large scale EM modes such as the aforementioned LLM. Finally in section 6, the mechanisms governing the ITB formation in different stages are concluded and a full picture of ITB’s self-sustainment is proposed.

## II. SIMULATION SETUP

All simulations reported in this paper are performed with the first-principle gyrokinetic code GENE [27] in flux-tube version. The simulated flux-tube is at $\rho_{tor,0} =$
<table><caption>TABLE I: Input parameters derived from HL-2A discharge #22453 at two experimental time 510 and 650 ms respectively. Parameters of electron are used for reference. For 510 ms, $n_e=1.272 \times 10^{19} m^{-3}$, $T_e=0.811$ keV, and $B_t = 1.349$ T, while for 650 ms, $n_e=1.271 \times 10^{19} m^{-3}$, $T_e=0.853$ keV, and $B_t = 1.345$ T. The aspect ratio R/a=4.87 and major radius R=1.68 m. In GENE code, electron $\beta_e = 8\pi n_e T_e / B_t^2$ serves as a reference value and the total $\beta = (1 + \sum_{j \neq e} n_j T_j / n_e T_e) \beta_e$. For the case without fast deuterium ion, its density $n_{fi}$ is set to zero and that of thermal deuterium ion, $n_i$, is set equal to $n_e$ to ensure quasi-neutrality.</caption><thead><tr><th>t<sub>exp</sub>[ms]</th><th>n<sub>i</sub>/n<sub>e</sub></th><th>n<sub>fi</sub>/n<sub>e</sub></th><th>T<sub>i</sub>/T<sub>e</sub></th><th>T<sub>fi</sub>/T<sub>e</sub></th><th>R/L<sub>n<sub>e</sub></sub></th><th>R/L<sub>n<sub>i</sub></sub></th></tr></thead><tbody><tr><td>510</td><td>0.75</td><td>0.25</td><td>0.84</td><td>22.62</td><td>4.36</td><td>5.68</td></tr><tr><td>650</td><td>0.72</td><td>0.28</td><td>1.29</td><td>22.22</td><td>6.83</td><td>9.13</td></tr></tbody><thead><tr><th>t<sub>exp</sub>[ms]</th><th>R/L<sub>n<sub>fi</sub></sub></th><th>R/L<sub>T<sub>e</sub></sub></th><th>R/L<sub>T<sub>i</sub></sub></th><th>R/L<sub>T<sub>fi</sub></sub></th><th>β<sub>e</sub>[%]</th><th>γ<sub>E</sub></th></tr></thead><tbody><tr><td>510</td><td>0.40</td><td>7.00</td><td>9.53</td><td>1.47</td><td>0.23</td><td>0.065</td></tr><tr><td>650</td><td>0.80</td><td>7.22</td><td>25.61</td><td>1.67</td><td>0.24</td><td>0.187</td></tr></tbody></table>

0.25, slightly inside the ITB foot. Here the subscript '0' indicates flux-tube location. Miller geometry [28] is extracted from EFIT equilibrium. An extended region of low but positive magnetic shear $\hat{s}$ is observed inside the ITB foot, and at the simulated location $\hat{s} = 0.12$ with the safety factor $q_0 = 1.05$. Typical grid parameters are as follows: perpendicular box sizes [$L_x, L_y$] = [272, 218] in units of ion Larmor radius $\rho_i$ with discretizations [$n_x, n_y$] = [768, 48], $n_z = 32$ points in parallel direction, 32 points in parallel velocity directions and 32 magnetic moments. Here, x is the radial coordinate defined as $x = a\rho_{tor}$ (a the minor radius), y the binormal coordinate and z the coordinate along the field line. When the effects of the perpendicular flow shear are considered, large aspect ratio and circular poloidal cross-section are assumed and therefore the normalized mean $E \times B$ shearing rate is defined as $\gamma_E \equiv (\rho_{tor,0}/q_0)(d\Omega/d\rho_{tor})/(c_s/R)$. $\Omega$ is the toroidal angular velocity, R the major radius and $c_s$ the sound speed. The full impact of $\gamma_E$ is considered in non-linear simulations only, in order to ensure the compatibility of the $E \times B$ algorithm [29] implemented in GENE. Other physical parameters are shown in Table I. With the aim of analyzing the individual effects of fast ion, finite-$\beta$ and $\gamma_E$, simulations are divided into subsets with or without some of these parameters, and they are performed at both the ITB triggering time, 510 ms, and when the ITB is well-developed at 650 ms. For the case without fast deuterium ion, its density $n_{fi}$ is set to zero and the thermal ion $n_i$ is set equal to $n_e$ to ensure quasi-neutrality. Note that the tuple $(n_{fi}, \beta, \gamma_E)$ is frequently used in the following figures to indicate the simulation conditions.

### III. INSTABILITIES

The frequencies and growth rates of the most unstable modes in linear simulations are presented in Fig. 2. At both 510 and 650 ms, the spectra are dominated by the electrostatic (ES) ITG modes, which are characterized

by frequencies in direction of ion diamagnetic drift and peaks at binormal wave number $k_y \rho_i \approx 0.3$ (or equally toroidal number $n \approx 12$). It can be observed, comparing the cases with and without $\beta$, that ITG is stabilized by the well-known linear finite-$\beta$ effects [30, 31]. Furthermore, fast ions exert another damping effect on ITG. This damping effect arises from the dilution of main ion species, which act as the driven force of ITG, and thus reduces ITG growth rates by a factor scaling with the fast-ion concentration $n_{fi}/n_e$ [11, 32]. To confirm this, $n_{fi}/n_e$ is scanned for the nominal case at 510 ms and the results in Fig. 3 clearly show that the thermal-ion-driven ITG, while insensitive to other characteristic parameters of fast ion, is stabilized by about 20% as $n_{fi}/n_e$ is increased from 0 to its nominal value. After the ITB is well developed at 650ms, modes at toroidal number n=1 and n=2, with frequencies higher than those of ITG modes, are found destabilized without the contribution of the fast ions but rather as a consequence of the combined influence of steep $R/L_{T_i}$, low $\hat{s}$ and finite-$\beta$. It was shown [33] previously that these EM modes have linear properties in good agreement with those of Beta-induced Alfvén Eigenmode (BAE) [34]. As analyzed by linear simulations therein [33], these BAEs are mainly destabilized by the thermal ion temperature gradient with the critical value as

$$R/L_{T_i}|_{\text{critical}} = \frac{1}{q_0 \sqrt{7/4 + T_e/T_i}} \frac{\omega_{tr}}{\omega_{*n_i}} \frac{R}{L_{n_i}}, \quad (1)$$

where $\omega_{tr} = \sqrt{2T_i/m_i/(q_0 R)}$ is the thermal ion transit frequency and $\omega_{*n_i}$ is the density part of the ion diamagnetic drift frequency $\omega_{*pi} = k_y T_i (R/L_{n_i} + R/L_{T_i})/(eB_t R)$. In brief, the distinct property of these modes is that their frequencies scale with both the transit and diamagnetic drift frequency of thermal ion, $\omega \sim \omega_{tr} \sim \omega_{*pi}$, nearly independent on the characteristic parameters of fast ion. Destabilized above a relatively low critical $\beta$, their mode structures in ballooning representation, unlike that of ITG modes which localized within small ballooning angle, have not only a large extension over the ballooning angle but also small scale variations with characteristic length of the order of $\beta^{1/2}$.

To gain an insight of the nonlinear characteristic of the instabilities, Fourier transforms are applied to fluctuating ES potential $\phi_1$ in the saturated phase of nonlinear simulations, and the results of several typical cases are averaged spatially and presented in Fig. 4. The zonal component of $\phi_1$, being the most prominent modes with zero frequency, are neglected in these spectra to highlight modes at $n \neq 0$. Two patterns of spectrum are generally observed comparing Fig. 4(a) and (b). For those cases where finite-$\beta$ and large $R/L_{T_i}$ are not jointly present, the nonlinear spectra are consistent with the linear results. As shown by the representative case in Fig. 3 (a), the peaks of the Fourier amplitudes coincide with the linear frequencies of most unstable modes, with bandwidths arising from nonlinear scattering of dominant modes or coexisting subdominant ones. For such
<figure><img src="image_3.png" /><figcaption>FIG. 2: Spectrum of growth rates (a) and frequencies (b) in linear simulations for the cases with ($n_{fi}$, $\beta$, $\gamma_E$) labeled above; here, 'n.' indicates the quantity is set to its nominal values in Table I and such notation will be used in the following figures.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3: Growth rates of the ITG at $n = 18$ (the peak for the nominal case at 510 ms) when fast ion density is scanned. Except for the fast ion density, other characteristic parameter of fast ion that has been changed in these scans is labelled in the legend. Note that density and density gradient of thermal ion have to be changed according to the quasi-neutrality condition.</figcaption></figure>

cases, no modes are found prominent in frequency range different from those of ITG modes, except that the bandwidth, as can be seen in Fig. 4(c) is broadened with finite $\gamma_E$. However, when finite-$\beta$ is considered in the presence of large $R/L_{T_i}$, high frequency modes appear apart from the ITGs. Those at $n=1$ and $n=2$ are the aforementioned BAEs with frequencies of around 40 and 55 kHz respectively, while those at higher $n$ have exponentially low amplitudes. It can be seen that, if the rotation frequency $f_{tor} \approx 6.7$ kHz is considered with $f = f_{lab} - n f_{tor}$, BAE at $n=2$ with $f \approx 55$ kHz corresponds to the mildly destabilized perturbation $f_{lab} \approx 70$ kHz in Fig. 1. When $\gamma_E$ is retained, as can be seen in Fig. 4(d) where the Fourier spectrum are averaged further over toroidal numbers, modes at $n=1$ appear with frequencies of $f \approx 14$ kHz, very close to the frequency of $n=1$ LLM in stationary frame. While finite $\gamma_E$ is essential for LLM to appear on one hand, fast ions act as a non-resonant energy source for its destabilization [15, 35] on the other. As can be seen in Fig. 4(d), $n=1$ LLM could appear but would be dominated by $n=1$ BAE, if the contribution of fast ions were excluded. Most importantly, the presence of a significantly large $R/L_{T_i}$ is indispensable for the destabilization of LLM.

## IV. ION HEAT FLUXES

The toroidal spectrum of the flux-surface-averaged ion heat fluxes, computed from the saturated phase of nonlinear simulations in unit of gyroBohm (gB $\equiv n_e T_e^{5/2} m_i^{1/2} / (e^2 B_t^2 a^2)$), are presented in Fig. 5(a) and (c). For comparison, they are integrated over toroidal number $n$ and plotted in the right panels as a function of the quasi-linear ion heat fluxes calculated from the corresponding linear growth rates according to [36]. In addition, a model constant $C$ is determined by the ratio of non-linear to quasi-linear flux level of the full case at 510 ms, and a dashed line indicating the prediction of quasi-linear model $Q_i = C Q_i^{qs}$ is shown in Fig. 5(b) and (d). Quasi-linear theory contains no information of the exact nonlinear couplings between instabilities, typically the couplings with zonal components. Therefore, comparison between quasi-linear and non-linear fluxes disentangles the non-linear behaviours of the instabilities from the linear ones, and thus serves as an indicator of the nonlinear effects. In the same panels, the ion heat fluxes calculated by the transport code ONETWO [37] are shown by dash-dotted lines to indicate the power balance level. Also, the heat fluxes of fast ions are typically negligible compared to those of thermal ions, and therefore omitted in the following analysis. At 510 ms, it can be seen from Fig 5.(a) that the fluxes are significantly reduced when fast ions and/or finite-$\beta$ are included and that the effect of fast ions is more effective for the reduction. Comparing the full case and the case without any factor, the fluxes are found to be reduced by about 78%, to which fast ion alone contribute about 90%. Fast ions in our cases have large density but relatively low pressure gradient, destabilizing no extra mode as analyzed above. Their introduction in our simulations merely cause a dilution of the fraction of thermal ions, the latter being the main drive of ITG responsible for most of the fluxes. Such fast-ion dilution effect is closely pertinent to tokamak with low plasmas density and high NBI power, and is recognized as the key factor for the triggering of ITB at 510 ms. Another noteworthy point is that, as can be seen in Fig. 5(b), the non-linear fluxes are highly pre-
<figure><img src="image_5.png" /><figcaption>FIG. 4: Frequency spectra of $\phi_1$ in nonlinear simulations. Logarithmic $\langle|F[\phi_1]|^2\rangle_{k_x,z}$ for the case with ($n_{fi}, \beta, \gamma_E$) = (n, n, 0.) at (a) 510 ms and (b) 650 ms; red dotted line indicating the linear frequencies are aslo shown for the respective case in (a) and (b);normalized $\langle|F[\phi_1]|^2\rangle_{k_x,k_y,z}$ at (c) 510 ms and (d) 650 ms for cases as labeled in the legend.</figcaption></figure>

dictable by quasi-linear theory from the growth rates in the corresponding linear cases. This indicates that the mechanisms behind the flux reduction at 510 ms mainly lies in the linear stabilization effects of both fast-ion di- lution and finite-$\beta$ as shown in Fig. 2 (a).

The quasi-linear prediction are reliable also at 650 ms before finite-$\beta$ is included. In Fig. 5(c), ITG modes are fully destabilized due to steep $R/L_{T_i}$ and the fluxes driven by them, peaking around n=6 or $k_y\rho_i = 0.15$, would reach 50 gB in total if all the stabilizing factors were discarded. The effect of fast-ion dilutions are not sufficient to suppress such fluxes. Instead, it is only when finite-$\beta$ is retained that the fluxes are drastically reduced. More importantly, the flux reduction caused by finite-$\beta$ in nonlinear simulations largely overtakes what is pre- dicted by quasi-linear model. It is seen in Fig. 4(d) that the effect of finite-$\beta$ alone can reduce the total 50 gB fluxes to 7.5 gB, deviating from the quasi-linear pre- diction by about 74%. When all the other factors are considered along with finite-$\beta$, the total fluxes are even- tually reduced to around 2.7 gB. Non-linear finite-$\beta$ ef- fect is thus identified as the dominant stabilizing effects during the sustainment phase of ITB. The underlying mechanism was investigated previously in [36] and at- tributed to an energy transferring enhanced by finite-$\beta$ between the flux-driven ITG modes and zonal flows. As will be shown below, a link between the flux reduction and the prominent growth of zonal flows is indeed identi- fied with finite-$\beta$, but the energy required for such growth are mainly tapped from n=1 EM modes instead of the ITG ones.

The inclusions of both fast ion and finite-$\beta$ have made a significant contribution to drawing the simulated heat fluxes near the power balance fluxes, but finite difference still exist between them at both 510 and 650 ms. Such differences may arise from the inevitable errors in mea- surement and parameters evaluations, but they could do little harm to our conclusions which rely mainly on the relative change of fluxes rather than on their absolute val- ues under different simulation conditions. By sharp con- trast to the beneficial role played by fast ion and finite-$\beta$, the effects of $E \times B$ shearing on the fluxes are barely vis- ible at both 510 and 650 ms. It is found in Fig. 5 that the retaining of finite $\gamma_E$ slightly increased the fluxes, but such change is within the error bar. The ineffective- ness of $\gamma_E$ was reported also in other tokamaks [38, 39] and can be simply explained by its low value, which in our case is only about one third of the growth rates of dominating ITG modes (if compared in the same unit $c_s/a$). Therefore, it is concluded that the $E \times B$ shear- ing has not direct impact on the ITB formation other than changing the characteristics of the n=1 EM modes, which, in spite of their dominant amplitudes in frequency spectrum, drive much less fluxes than their ITG compan- ions.

V. FLUX REDUCTION AND ZONAL FLOWS

The aforementioned discrepancies between non-linear and quasi-linear fluxes with finite-$\beta$ are related to the effect of the zonal flows (ZF). To confirm this, the full
<figure><img src="image_6.png" /><figcaption>FIG. 5: Spectra of ion heat fluxes for the cases at time 510 (a) and 650 ms (c) with $(n_{fi}, \beta, \gamma_E)$ as labeled above; the spectra at $n > 24$ are omitted for they are low in values; total non-linear fluxes v.s. quasi-linear fluxes at time 510 (b) and 650 ms (d). The levels of ion power balance heat flux calculated by ONETWO, 0.58 and 0.54 gB, are also shown in (b) and (d) respectively.</figcaption></figure>

time traces of total ion fluxes are shown in Fig. 6(a) for the cases with and without finite-$\beta$ at 650 ms, labelled by EM and ES respectively, and in the same plot the ZF energies for the respective case are also displayed. Here, the field energy of each $n$ is defined concerning only the ES part as $E_n \equiv \sum_{k_x} \int JdzC_1|\phi_1|^2$, where $J$ is the Jacobian. A positive-defined real constant $C_1(k_x, k_y, z)$ is included so that $E_n$ corresponds to the field part of the free energy [40], and could be substituted with other positive-defined real constant such as $k_\perp^2$. From Fig. 6(a), it's observed regardless of whether finite-$\beta$ is retained or not, that the ion heat fluxes at first develop linearly to form the $\gamma_n$-dependent shape during the initial phase, the time window before $t_1$ when the amplitudes of ZFs are low, and that the peaks of the spectra begin to shift toward lower toroidal number as the ZFs continue to develop in the transient phase from $t_1$ to $t_2$. The corresponding flux spectra are shown in Fig. 6(b). At this phase, the heat flux spectra tend to evolve into the $\gamma_n/k_\perp^2$-dependent [41, 42] quasi-linear shape with overall values predictable from the dashed line in Fig 5(d). For the ES case, the heat fluxes begin to saturate at the quasi-linear level. For the EM case, however, it is observed that the heat fluxes, instead of becoming saturated at the quasi-linear level, continue to abate slowly as the ZF energy in such case is experiencing a persistent growth. Since the quasi-linear model predicts ion heat flux without any information of ZF, the effect of ZF growth can be simply estimated by the deviation between the saturated nonlinear heat flux and the quasi-linear one (as indicated in Fig. 5(d)), while the cause of such ZF growth is investigated in the following. The time evolution of $E_n$ naturally depends

on the linear contributions and nonlinear ones, but only the latter contribute to the net growth of ZF energy $E_0$. Taking the time derivative of $E_n$ and substituting $\phi_1$ with the modified distribution function $g_1$ through the field equation, the nonlinear contribution to $dE_n/dt$ is expressed as (see the Appendix)

$$ \frac{dE_n}{dt}|_{NL} \equiv \sum_{n'} T(n, n', t) \\ = Re \sum_{k'_y, k_x, k'_x} (k'_x k_y - k'_y k_x) \int Jdz M \bar{\phi}_1^* |k \chi_{1j} | k' g_{1j} | k-k', \quad (2) $$

where $M = \sum_j \pi n_j q_j \int dv_\parallel d\mu B_0$ is a moment operator and $\chi_1 = \bar{\phi}_1 - v_{th,j} v_\parallel \bar{A}_{1,\parallel} + T_j \mu / q_j \bar{B}_{1,\parallel}$ the gyro-averaged effective potential. $T(n, n', t)$ is calculated focusing on the coupling between ZF ($n=0$) and all the other $n'$, and the results shown in Fig. 6(c) have been averaged over the time window when there is a net growth of ZF. It is thus seen that ZF has mainly drained energies from $n=8 \sim 12$ ITG components when the low-n EM modes are artificially suppressed by neglecting finite-$\beta$ effect. Instead, when finite-$\beta$ is retained and the low-n EM modes are destabilized by the large $R/L_{T_i}$, ZF receives a significant positive portion of energy from these low-n modes, mostly from $n=1$ LLM, and develops to a much larger amplitude than that in the case without finite-$\beta$. Such favorable energy transfer is an evidence of the self-regulatory system where an EM mode, serving as a catalyst, transfers the free energy it obtained from the ITB-generated large $R/L_{T_i}$ to the ZF which helps mitigate the heat fluxes and in turn sustain the ITB. Consequently, a self-organized mechanism is proposed which
is characterized by an energy transfer that is facilitated by the saturation of the low-n EM modes, in this case the LLM, and that results in the increase of ZF activities and reduction of heat fluxes.

## VI. DISCUSSION AND CONCLUSIONS

The ITB characteristics in HL-2A have been analyzed by performing non-linear gyrokinetic simulations. The emphasis of our study have been placed on the effects of fast ions, finite-β and $E \times B$ shear. It is found that the complete ITB formation process can be conceptually divided into two stages where distinct mechanisms dominate. Widely effective as it is, the $E \times B$ shear stabilization in our cases is not found to play a remarkable role in any of these stages, mainly because the shearing rate $\gamma_E$ is relatively low compared to the ITG growth rates. In the first stage, the plasmas instabilities are dominated by ITG modes which are subjected to the stabilizing effects of both finite-β and fast ions. It is found that the triggering of ITB is mainly caused by the stabilization of ITG under the effect of fast-ion dilution which basically depend on linear physics. Once the ITB is fully developed, the sustaining of the ITB is determined by the reduction of heat turbulent transport by large scale zonal flows. In this second stage, the steep ITB-generated pressure gradient, combined with the effect of finite-β, is able to bring about an abundant varieties of large scale EM modes, in our case the LLM. Instead of driving significant fluxes, LLM acts as a catalyst that transfers the ITB free energy obtained during the triggering process to the zonal flows, which in turn mitigate the flux and sustain the ITB ultimately. The full ITB formation is therefore characterized as a self-regulatory multi-scale physics system leading to a self-sustained ITB. Although these conclusions are obtained from simulations which employ several simplifying model, such as the local assumption and Maxwellian fast ions, they provide an initial picture of the process of ITB formations. It has been shown [43] that the stabilization effects of NBI fast ions, such as the dilution, are not significantly affected by their distribution in the velocity space. But to validate our conclusions in further, more works are needed employing realistic fast ion distribution, and the global gyrokinetic simulations that is much more demanding computationally may be necessary to rigorously account for the effect of large scale flow shear and to completely accommodate all the modes involved. Besides, similar study should be carried out to HL-2A plasmas with FB and ITB, as the correlation found here between LLM and ITB may also exist between FB and ITB. Nevertheless, the mechanism proposed in this paper could be important if LLM can be induced, e.g. by tailoring the q-profile, to future tokamak devices like ITER with low $E \times B$ shearing, which is not found to play a major role here on any stage of ITB formation.

## VII. ACKNOWLEDGEMENT

The authors are very grateful to Mr. Chen Qian, Mr. Zhang Xing, Mr. Fang Kairui, Dr. Hao Guangzhou, Dr. Yu Deliang and the HL-2A experiment team for providing and processing experimental data. This work was supported by the National Natural Science Foundation of China with grant Nos. 12275071 and U1967206 and also partially by National Key R&D Program of China under Grant Nos. 2017YFE0301200 and 2017YFE0301201.

## VIII. APPENDIX

The derivations of Eq. 2 are reported in the following. In flux-tube version of GENE, the normalized gyrokinetic Vlasov equation for the modified distribution function $g_{1j}$ of species j can be written as

$$ \frac{dg_{1j}}{dt} = L_G\chi_{1j} + L_C(g_{1j} + q_j\chi_{1j}\frac{F_{0j}}{T_j}) + L_\parallel(g_{1j} + q_j\chi_{1j}\frac{F_{0j}}{T_j}) - \{\chi_{1j}, g_{1j}\}_{x,y}, \quad (3) $$

where $L_G$ is the gradient prefactor,

$$ L_G = \left( \left( -\frac{3}{2} + v_\parallel^2 + \mu B_0 \right) \frac{R}{L_{T_j}} + \frac{R}{L_{n_j}} \right) F_{0j} ik_y, \quad (4) $$

$L_C$ the curvature prefactor,

$$ L_C = \frac{T_j(2v_\parallel^2 + \mu B_0)}{q_j B_0} K_x ik_x + \left( -\frac{T_j(2v_\parallel^2 + \mu B_0)}{q_j B_0} K_y + \beta \frac{T_j v_\parallel^2 p_0}{q_j B_0^2} \frac{R}{L_{p_0}} \right) ik_y, \quad (5) $$

and $L_\parallel$ the parallel-dynamic operator,

$$ L_\parallel = v_{th,j} \frac{F_{0j}}{2} \frac{1}{JB} \left\{ \frac{1}{F_{0j}}, \right\}_{v_\parallel,z}. \quad (6) $$

Here, the Poisson bracket of two arbitrary function f and g over the variables u and v is defined as

$$ \{f,g\}_{u,v} = \frac{\partial f}{\partial u} \frac{\partial g}{\partial v} - \frac{\partial f}{\partial v} \frac{\partial g}{\partial u}. \quad (7) $$

When the nonlinear term in Eq. 3 (the last term in the right hand side) is evaluated in Fourier space at $(k_x, k_y, z)$, the multiplications in the Poisson bracket are transformed into convolutions, i.e.

$$ \begin{aligned} & - \{\chi_{1j}, g_{1j}\}_{x,y|k} \\ &= \sum_{k'_x, k'_y} (k'_x k_y - k_x k'_y) \chi_{1j}|_{k'|g_{1j}|_{k-k'}}. \end{aligned} \quad (8) $$

Full magnetic fluctuations are considered in $\chi_1 = \bar{\phi}_1 - v_{th,j}v_\parallel \bar{A}_{1,\parallel} + T_j\mu/q_j \bar{B}_{1,\parallel}$, the gyro-averaged effective potential, where the bar over quantities indicates gyro-average. Note that in local limit the gyro-average of ES
<figure><img src="image_7.png" /><figcaption>FIG. 6: (a) Time trace of the total ion fluxes and ZF energy, labelled by $Q_i$ and ZF respectively, for the cases at 650 ms with ($n_{fi}$, $\beta$, $\gamma_E$) equal to (n, 0, n.) and (n, n, n.), labelled by EM and ES respectively. (b) Flux spectrum at the normalized time marked in (a). (c) Time-averaged nonlinear term contribution from each $n'$ to the growth of ZF energy for the EM and ES case during the labelled normalized time.</figcaption></figure>

potential $\phi_1$ is simply $\bar{\phi}_1 = J_0(k_\perp \rho_j) \phi_1$, where the $J_n$ is the Bessel function of $n^{th}$ order. The symbol $|_k$ indicates the quantity before it is evaluated at ($k_x, k_y$). In the curvature term, the $K_x$ and $K_y$ are the curvature factor in radial and binormal direction respectively. Their definition, as well as those of other quantities, can be found in, e.g. ref. [44] and [45]. The field equation of the ES potential $\phi_1$ is coupled with that of the parallel fluctuating magnetic field $B_{1\parallel}$ when finite-$\beta$ is considered. The coupled field equations are

$$C_1 \phi_1 + C_2 B_{1\parallel} = M J_0(k_\perp \rho_j) g_{1j}, \quad (9)$$

$$C_2 \phi_1 + C_3 B_{1\parallel} = M \frac{2J_1(k_\perp \rho_j) T_j \mu}{k_\perp \rho_j} g_{1j}, \quad (10)$$

from which we obtain

$$\phi_1 = \frac{C_3}{C_1 C_3 - C_2^2} M J_0(k_\perp \rho_j) g_{1j} - \frac{C_2}{C_1 C_3 - C_2^2} M \frac{2J_1(k_\perp \rho_j) T_j \mu}{k_\perp \rho_j} g_{1j}. \quad (11)$$

The moment operator $M$ is defined as

$$M = \sum_j \pi n_j q_j \int dv_{\parallel} d\mu B_0, \quad (12)$$

and the definitions of the coefficients $C_1$, $C_2$ and $C_3$ (which are real and only depend on $k_x$, $k_y$ and $z$) can be found in Page 33 of ref. [44]. By Eq. 11, the total derivative of the ES energy at $k_y$ can be expressed as

$$\begin{aligned} \frac{dE_n}{dt} &= \sum_{k_x} \int J dz C_1 \left( \frac{\partial \phi_1^*}{\partial t} \phi_1 + \phi_1^* \frac{\partial \phi_1}{\partial t} \right) \\ &= Re \sum_{k_x} \int J dz C_1 \phi_1^* \frac{\partial \phi_1}{\partial t} \\ &= Re \left( \sum_{k_x} \int J dz \frac{C_1}{C_1 C_3 - C_2^2} \phi_1^* \right. \\ & \qquad \left. \times M(C_3 J_0(k_\perp \rho_j) - C_2 \frac{2J_1(k_\perp \rho_j) T_j \mu}{k_\perp \rho_j}) \frac{\partial g_{1j}}{\partial t} \right). \end{aligned} \quad (13)$$

As the nonlinear contribution alone is of our concern, we can obtain, by substituting the $\partial g_{1j}/\partial t$ in Eq. 13 with only the nonlinear term Eq. 8,

$$\begin{aligned} \frac{dE_n}{dt}|_{NL} = & Re \sum_{k'_y, k_x, k'_x} (k'_x k_y - k'_y k_x) \int J dz M \phi_1^*|_k \\ & \times \left( \frac{C_1 C_3}{C_1 C_3 - C_2^2} J_0(k_\perp \rho_j) \right. \\ & \left. - \frac{C_1 C_2}{C_1 C_3 - C_2^2} \frac{2J_1(k_\perp \rho_j)}{k_\perp \rho_j} \right) \\ & \times \chi_{1j} |k' g_{1j}|_{k-k'}, \end{aligned} \quad (14)$$

where the commutativity between the moment operator $M$ and any spatial quantities has been used. As the coefficient $C_3$ is inversely proportional to $\beta$ which is closed to zero for our case, the second term in Eq. 14 makes negligible contribution to the total value. Therefore, taking
the limit $C_3 \to \infty$, Eq. 2 can be obtained from Eq. 14.

But note that Eq. 14 instead of 2 is actually used to
produce Fig. 6 (c) for the sake of completeness.

[1] F. Romanelli, Physics of Fluids B: Plasma Physics 1, 1018 (1989).

[2] J. W. Connor, T. Fukuda, X. Garbet, C. Gormezano, V. Mukhovatov, M. Wakatani, the ITB Database Groupa, and the ITPA Topical Group on Transport and Internal Barrier Physics, Nucl. Fusion 44, R1 (2004).

[3] R. E. Bell, F. M. Levinton, S. H. Batha, E. J. Synakowski, and M. C. Zarnstorff, Plasma Phys. Controlled Fusion 40, 609 (1998).

[4] M. Romanelli, A. Zocco, and F. Crisanti, Plasma Phys. Controlled Fusion 52, 045007 (2010).

[5] E. J. Doyle, C. M. Greenfield, M. E. Austin, L. R. Baylor, K. H. Burrell, T. A. Casper, J. C. DeBoo, D. R. Ernst, C. Fenzi, P. Gohil, R. J. Groebner, W. W. Heidbrink, G. L. Jackson, T. C. Jernigan, J. E. Kinsey, L. L. Lao, M. Makowski, and et al., Nucl. Fusion 42, 333 (2002).

[6] E. Joffrin, C. D. Challis, G. D. Conway, X. Garbet, A. Gude, S. Günter, N. C. Hawkes, T. C. Hender, D. F. Howell, G. T. A. Huysmans, E. Lazzaro, P. Maget, M. Marachek, A. G. Peeters, S. D. Pinches, S. E. Shara-pov, and JET-EFDA contributors, Nucl. Fusion 43, 1167 (2003).

[7] Y. Sakamoto, T. Suzuki, S. Ide, Y. Koide, H. Takenaga, Y. Kamada, T. Fujita1, T. Fukuda, T. Takizuka, H. Shiraï, N. Oyama, Y. Miura, the JT-60 Team, K. W. Hill, and G. Rewoldt, Nucl. Fusion 44, 876 (2004).

[8] D. L. Yu, Y. L. Wei, L. Liu, J. Q. Dong, K. Ida, K. Itoh, A. P. Sun, J. Y. Cao, Z. B. Shi, Z. X. Wang, Y. Xiao, B. S. Yuan, H. R. Du, X. X. He, W. J. Chen, Q. Ma, and et al., Nucl. Fusion 56, 056003 (2016).

[9] P. H. Diamond, V. B. Lebedev, D. E. Newman, B. A. Carreras, T. S. Hahm, W. M. Tang, G. Rewoldt, and K. Avinash, Phys. Rev. Lett. 78, 1472 (1997).

[10] K. H. Burrell, Phys. Plasmas 4, 1499 (1997).

[11] G. Tardini, J. Hobirk, V. Igochine, C. Maggi, P. Martin, D. McCune, A. Peeters, A. Sips, A. Stäbler, J. Stober, and the ASDEX Upgrade Team, Nucl. Fusion 47, 280 (2007).

[12] A. D. Siena, R. Bilato, T. Görler, A. Bañon-Navarro, E. Poli, V. Bobkov, D. Jarema, E. Fable, C. Angioni, Y. O. Kazakov, R. Ochoukov, P. Schneider, M. Weiland, F. Jenko, and the ASDEX Upgrade Team, Phys. Rev. Lett. 127, 025002 (2021).

[13] Y. L. Wei, D. L. Yu, L. Liu, K. Ida, M. v. Hellermann, J. Y. Cao, A. P. Sun, Q. Ma, W. J. Chen, Y. Liu, L. W. Yan, Q. W. Yang, X. R. Duan, and Y. Liu, Review of Scientific Instruments 85, 103503 (2014).

[14] P. H. Diamond, S. I. Itoh, K. Itoh, and T. S. Hahm, Plasma Phys. Controlled Fusion 47, R35 (2005).

[15] R. B. Zhang, X.-Q. Wang, C.-J. Xiao, X.-G. Wang, Y. Liu, W. Deng, W. Chen, X.-T. Ding, X.-R. Duan, and the HL-2 A Team, Plasma Phys. Controlled Fusion 56, 095007 (2014).

[16] L. Yu, W. Chen, X. Ding, X. Ji, Z. Shi, D. Yu, M. Jiang, D. Li, J. Li, Y. Li, Y. Zhou, R. Ma, W. Li, B. Feng, Y. Huang, X. Song, and et al., Journal of the Physical

Society of Japan 86, 024501 (2017).

[17] S. Günter, A. Gude, J. Hobirk, M. Maraschek1, S. Saarelma, S. Schade, R. C. Wolf, and A. U. Team, Nucl. Fusion 41, 1283 (2002).

[18] W. Ge, Z.-X. Wang, F. Wang, Z. Liu, and L. Xu, Nucl. Fusion 63, 016007 (2023).

[19] S. Wang, H. Cai, X. Chen, and D. Li, Plasma Phys. Controlled Fusion 65, 055018 (2023).

[20] X. X. He, L. W. Yan, D. L. Yu, W. Chen, L. M. Yu, Q. Ma, L. Liu, Y. L. Wei, X. F. He, N. Zhang, Y. G. Li, T. F. Sun, L. Nie, Z. B. Shi, W. L. Zhong, M. Xu, and X. R. Duan, Plasma Phys. Controlled Fusion 64, 015007 (2022).

[21] W. Deng, Y. Liu, W. L. Ge, M. Jiang, Z. B. Shi, D. Li, X. Q. Ji, Y. B. Dong, F. Wang, J. Y. Cao, W. L. Zhong, J. M. Gao, Y. P. Zhang, G. Y. Li, Z. X. Wang, M. Xu, and X. R. Duan, Phys. Plasmas 29, 102106 (2022).

[22] L. M. Yu, W. Chen, M. Jiang, Z. B. Shi, X. Q. Ji, X. T. Ding, Y. G. Li, R. R. Ma, P. W. Shi, S. D. Song, B. S. Yuan, Y. Zhou, R. Ma, X. M. Song, J. Q. Dong, M. Xu, Y. Liu, L. W. Yan, Q. W. Yang, Y. H. Xu, D. X. R., and HL-2A Team, Nucl. Fusion 57, 036023 (2017).

[23] E. Joffrin, C. D. Challis, T. C. Hender, D. F. Howell, and G. T. A. Huysmans, Nucl. Fusion 42, 235-242 (2017).

[24] E. Joffrin, R. Wolf, B. Alper, Y. Baranov, C. D. Challis, M. d. Baar, C. Giroud, C. W. Gowers, N. C. Hawkes, T. C. Hender, M. Marachek, D. Mazon, V. Parail, A. Peeters, K.-D. Zastrow, and contributors to the EFDA-JET Workprogramme, Plasma Phys. Controlled Fusion 44, 1203-1214 (2002).

[25] E. Joffrin, G. Gorini, C. D. Challis, N. C. Hawkes, T. C. Hender, D. F. Howell, P. Maget, P. Mantica, D. Mazon, S. E. Sharapov, G. Tresset, and contributors to the EFDA-JET Workprogramme, Plasma Phys. Controlled Fusion 44, 1739-1752 (2002).

[26] K. Ida and T. Fujita, Plasma Phys. Controlled Fusion 60, 033001 (2018).

[27] F. Jenko, W. Dorland, M. Kotschenreuther, and B. N. Rogers, Phys. Plasmas 7, 1904 (2000).

[28] R. L. Miller, M. S. Chu, J. M. Greene, Lin-Liu, Y. R., and R. E. Waltz, Phys. Plasmas 5, 973 (1998).

[29] B. F. McMillan, J. Ball, and S. Brunner, Plasma Phys. Controlled Fusion 61, 055006 (2019).

[30] J. Weiland and A. Hirose, Nucl. Fusion 32, 151 (1992).

[31] A. Hirose, Phys. Plasmas 7, 433 (2000).

[32] G. J. Wilkie, A. Iantchenko, I. G. Abel, E. Highcock, I. Pusztai, and J. Contributors, Nucl. Fusion 58, 082024 (2018).

[33] W. H. Lin, J. Q. Li, J. Garcia, and S. Mazzi, Chinese Physic B 32, 025202 (2022).

[34] F. Zonca, L. Chen, and R. Santoro, Plasma Phys. Controlled Fusion 38, 2011 (1996).

[35] W. Xie, Physical Review E 105, 055208 (2022).

[36] G. Whelan, M. Pueschel, and P. Terry, Phys. Rev. Lett. 120, 175002 (2018).

[37] F. B. Pfeiffer, W. and Marcus, C. J. Armentrout, G. L. Jahns, T. W. Petrie, and R. E. Stockdale, Nucl. Fusion
25, 655 (1985).

[38] J. Citrin, J. Garcia, T. Görler, F. Jenko, P. Mantica, D. Told, C. Bourdelle, D. Hatch, G. Hogeweij, T. Johnson, et al., Plasma Physics and Controlled Fusion 57, 014032 (2014).

[39] C. Pan, G. M. Staebler, L. L. Lao, A. M. Garofalo, X. Gong, Q. Ren, J. McClenaghan, G. Li, S. Ding, J. Qian, B. Wan, G. S. Xu, W. Solomon, M. O., and S. P. Smith, Nucl. Fusion 57, 036018 (2017).

[40] A. Bañon Navarro, P. Morel, M. Albrecht-Marc, D. Carati, F. Merz, T. Görler, and F. Jenko, Phys. Plasmas 18, 092303 (2011).

[41] T. D. F. Jenko and C. Angioni, Plasma Phys. Controlled Fusion 47, B195-B206 (2005).

[42] C. Bourdelle, X. Garbet, F. Imbeaux, A. Casati, N. Dubuit, R. Guirlet, and T. Parisot, Phys. Plasmas 14, 112501 (2007).

[43] A. Di Siena, T. Görler, H. Doerk, R. Bilato, J. Citrin, T. Johnson, M. Schneider, E. Poli, and J. Contributors, Physics of Plasmas 25, 042304 (2018).

[44] F. Merz, Gyrokinetic simulation of multimode plasma turbulence, Ph.D. thesis, Universität Münster (2008).

[45] D. Told, Gyrokinetic Microturbulence in Transport Barriers, Ph.D. thesis, University of Ulm (2012).