A self-sustaining mechanism for internal transport barrier formation in HL-2A tokamak plasmas

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 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 (EM) 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 EM 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][4][5][6][7][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 × 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 × 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 × 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 × 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 selfsustainment 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 ω = 0), which tap energy non-linearly from large scale MagnetoHydroDynamics (MHD) fluctuations that are destabilized by the ITB itself.Meanwhile, the E × 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=510ms, 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 Ti (≈ 20) with the ITB foot located at ρ tor ≈ 0.4.Here, R is the major radius, L Ti the inverse logarithmic gradient of ion temperature and ρ 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][18][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.The dynamic interplay between LLM and ITB remains obscure so far and will be investigated further in this work.The structure of this paper is arranged as follow: after the simulation setup is addressed in section 2, the dominant instabilities in various simulation conditions are an-alyzed 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 benefitted 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 on 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 [22] in flux-tube version.The simulated flux-tube is at ρ tor,0 = 0.25, slightly inside the ITB foot.Here the subscript '0' indicates flux-tube location.Miller geometry [23] is extracted from EFIT equilibrium.An extended region of low but positive magnetic shear ŝ is observed inside the Parameters of electron are used for reference.
For 510 ms, ne=1.272×10 19m −3 , Te=0.811 keV, and Bt = 1.349T, while for 650 ms, ne=1.271×10 19m −3 , Te=0.853 keV, and Bt = 1.345 T. The aspect ratio R/a=4.87 and major radius R=1.68 m.In GENE code, electron βe = 8πneTe/B 2 t serves as a reference value and the total β = (1 + j̸ =e n j T j /neTe)βe.For the case without fast deuterium ion, its density n f i is set to zero and that of thermal deuterium ion, n i , is set equal to ne to ensure quasi-neutrality.Here, x is the radial coordinate defined as x = aρ 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 ×B shearing rate is defined as γ E ≡ (ρ tor,0 /q 0 )(dΩ/dρ tor )/(c s /R).Ω is the toroidal angular velocity, R the major radius and c s the sound speed.The full impact of γ E is considered in non-linear simulations only, in order to ensure the compatibility of the E × B algorithm [24] implemented in GENE.Other physical parameters are shown in Table I.With the aim of analyzing the individual effects of fast ion, finite-β and γ 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.Note that the tuple (n f i , β, γ 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 ρ i ≈ 0.3 (or equally toroidal number n ≈ 12).It can be observed, comparing the cases with and without β, that ITG is stablized by the well-known linear finite-β effects [25,26].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 f i /n e [11,27].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 Ti , low ŝ and finite-β.It was shown [28] previously that these EM modes have linear properties in good agreement with those of Beta-induced Alfvén Eigenmode (BAE) [29].As analyzed by linear simulations therein [28], these BAEs are mainly destabilized by the thermal ion temperature gradient with the critical value as where ω tr = 2T i /m i /(q 0 R) is the thermal ion transit frequency and ω * ni is the density part of the ion diamagnetic drift frequency ω * pi = k y T i (R/L ni + R/L Ti )/(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, ω ∼ ω tr ∼ ω * pi , nearly independent on the characteristic parameters of fast ion.Destabilized above a relatively low critical β, 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 β 1/2 .
To gain an insight of the nonlinear characteristic of the instabilities, Fourier transforms are applied to fluctuating ES potential ϕ 1 in the saturated phase of nonlinear simulations, and the results of several typical cases are averaged spatially and presented in Fig. 3.The zonal component of ϕ 1 , being the most prominent modes with zero frequency, are neglected in these spectra to highlight modes at n ̸ = 0. Two patterns of spectrum are generally observed comparing Fig. 3(a) and (b).For those cases where finite-β and large R/L Ti 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 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. 3(c) is broadened with fintie γ E .However, when finite-β is considered in the presence of large R/L Ti , 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 ≈ 6.7 kHz is considered with f = f lab −nf tor , BAE at n=2 with f ≈ 55 kHz corresponds to the mildly destabilized perturbation f lab ≈ 70 kHz in Fig. 1.When γ E is retained, as can be seen in Fig. 3(d) where the Fourier spectrum are averaged further over toroidal numbers, modes at n=1 appear with frequencies of f ≈ 14 kHz, very close to the frequency of n=1 LLM in stationary frame.While finite γ E is essential for LLM to appear on one hand, fast ions act as a non-resonant energy source for its destabilization [15,30] on the other.
As can be seen in Fig. 3(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 Ti 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 t a 2 )), are presented in Fig. 4(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 [31].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 = CQ qs i is shown in Fig. 4(b) and (d).Quasi-linear theory contains no information of the exact nonlinear couplings between instabilities, typically the couplings with zonal components.Therefore, com- parison 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 [32] are shown by dash-dotted lines to indicate the power bal-ance 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 4 .(a)that the fluxes are significantly reduced when fast ions and/or finite-β are included and that the effect of fast ions is more effective for the reduc-tion.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. 4(b), the non-linear fluxes are highly predictable 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 dilution and finite-β as shown in Fig. 2 (a).The quasi-linear prediction are reliable also at 650 ms before finite-β is included.In Fig. 4(c), ITG modes are fully destabilized due to steep R/L Ti and the fluxes driven by them, peaking around n=6 or k y ρ 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-β is retained that the fluxes are drastically reduced.More importantly, the flux reduction caused by finite-β in nonlinear simulations largely overtakes what is predicted by quasi-linear model.It is seen in Fig. 4(d) that the effect of finite-β alone can reduce the total 50 gB fluxes to 7.5 gB, deviating from the quasi-linear prediction by about 74%.When all the other factors are considered along with finite-β, the total fluxes are eventually reduced to around 2.7 gB. Non-linear finite-β effect is thus identified as the dominant stabilizing effects during the sustainment phase of ITB.The underlying mechanism was investigated previously in [31] and attributed to an energy transferring enhanced by finite-β 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 identified with finite-β, 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-β 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 measurement 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 values under different simulation conditions.By sharp contrast to the beneficial role played by fast ion and finite-β, the effects of E × B shearing on the fluxes are barely visible at both 510 and 650 ms.It is found in Fig. 4 that the retaining of finite γ E slightly increased the fluxes, but such change is within the error bar.The ineffective-ness of γ E was reported also in other tokamaks [33,34] 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 × B shearing 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 companions.

V. FLUX REDUCTION AND ZONAL FLOWS
The aforementioned discrepancies between non-linear and quasi-linear fluxes with finite-β are related to the effect of the zonal flows (ZF).To confirm this, the full time traces of total ion fluxes are shown in Fig. 5(a) for the cases with and without finite-β 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 ≡ kx JdzC 1 |ϕ 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 [35], and could be substituted with other positive-defined real constant such as k 2 ⊥ .From Fig. 5(a) , it's observed regardless of whether finite-β is retained or not, that the ion heat fluxes at first develop linearly to form the γ 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. 5 (b).At this phase, the heat flux spectra tend to evolve into the γ n /k 2 ⊥ -dependent [36,37] quasi-linear shape with overall values predictable from the dashed line in Fig 4(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, continue to abate slowly as the ZF energy in such case is experiencing a persistent growth, whose cause is reported 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 ϕ 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) where M = j πn j q j dv ∥ dµB 0 is a moment operator and χ 1 = φ1 − v th,j v ∥ Ā1,∥ + T j µ/q j B1∥ 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. 5(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 ∼ 12 ITG components when the low-n EM modes are artificially suppressed by neglecting finite-β effect.Instead, when finite-β is retained and the low-n EM modes are destabilized by the large R/L Ti , 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-β.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/ Ti 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 × 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 × B shear stabilization in our cases is not found to play a remarkable role on any of these stages, mainly because the shearing rate γ E is ralatively low compared to the ITG growth rates.On 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.On 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.To validate our conclusions in further, 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.Nevertheless, the mechanism proposed in this paper could be important, e.g. if LLM can be induced, e.g. by tailoring the q-profile, to future tokamak devices like ITER with low E × 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 where L G is the gradient prefactor, L C the curvature prefactor, and L ∥ the parallel-dynamic operator, Here, the Poisson bracket of two arbitrary function f and g over the variables u and v is defined as 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.
Full magnetic fluctuations are considered in χ 1 = φ1 − v th,j v ∥ Ā1,∥ + T j µ/q j B1∥ , the gyro-averaged effective potential, where the bar over quantities indicates gyroaverage.Note that in local limit the gyro-average of ES potential ϕ 1 is simply φ1 = J 0 (k ⊥ ρ j )ϕ 1 , where the J n is the Bessel function of n th order.The symbol | k indicates the quatity 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. [38] and [39].The field equation of the ES potential ϕ 1 is coupled with that of the parallel fluctuating magnetic field B 1∥ when finite-β is considered.The coupled field equations are from which we obtain T j µ q j g 1j .
The moment operator M is defined as M = j πn j q j dv ∥ dµB 0 , 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. [38].By Eq. 11, the total derivative of the ES energy at k y can be expressed as As the nonlinear contribution alone is of our concern, we can obtain, by substituting the ∂g 1j /∂t in Eq. 13 with only the nonlinear term Eq. 8, where the commutativity between the moment operator M and any spatial quantities has been used.As the coefficient C 3 is inversely proportional to β 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 → ∞, Eq. 2 can be obtained from Eq. 14.
But note that Eq. 14 instead of 2 is actually used to produce Fig. 5 (c) for the sake of completeness.

FIG. 1 :
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 ρtor = 0.25 in the same discharge.

FIG. 2 :FIG. 3 :
FIG.2: Spectrum of growth rates (a) and frequencies (b) in linear simulations for the cases with (n f i , β, γ E ) labeled above; here, 'n.' indicates the quantity is set to its nominal values in TableIand such notation will be used in the following figures.

FIG. 4 :
FIG. 4: Spectra of ion heat fluxes for the cases at time 510 (a) and 650 ms (c) with (n f i , β, γ 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.

TABLE I :
Input parameters derived from HL-2A discharge #22453 at two experimetal time 510 and 650 ms respectively.