Global E × B flow pattern formation and saturation

The E × B flow staircase-like pattern observed in the first principle gyrokinetic numerical experiments of tokamak fusion plasmas forms due to a nonlinear time delay. Simulations demonstrate a finite time delay between the staircase occurrence in particle transport and that in the density profile. This novel finding shows that instability can arise from perturbations in transport and then influence the background turbulence. E × B flow staircase plays roles not only in shearing the transport but also as a nonlinear saturation mechanism of staircase instability. Experimental measurements in KSTAR tokamak L-mode plasmas are consistent with the numerical findings.


Introduction
Understanding the global pattern formation is a prominent subject in various research fields, such as potential vorticity (PV) staircases ( jets) in planetary atmospheres [1,2], thermohaline staircases in geophysics [3], zonal flow formation in the Earth's core [4], jams in traffic transportation [5] etc. In tokamak fusion plasmas, quasi-regular staircase-like corrugations are observed in the mean E × B (zonal) flow, strongly correlated pressure profiles and transport avalanches in numerical simulations [6,7] and experiments [8][9][10]. The mean E × B flow of interest here is the zonal (axisymmetric, fluxsurface averaged) component of E × B flows in a tokamak. Analogous to the strong influence of staircases in the general circulation in astrophysics and geophysics, as well as jams in transportation, the zonal flow staircase plays significant roles in the regulation of transport avalanches, and therefore is favorable for the enhancement of energy confinement in * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. tokamak fusion plasmas. Zonal flow staircases constrain the propagation of transport avalanches by generating localized, permeable transport barriers with enhanced local shearing layers through corrugating the global E × B pattern. This is a prime example of mesoscopic phenomena [11] with both academic and practical relevance.
E × B staircase-like corrugations are of great interest as they can favor the achievement of fusion energy. However, the underlying physical generation mechanism remains as conundrums to date. The heuristic model [12] proposed the staircase formation as being like a heat flux 'jam' that causes temperature profile corrugations, which resembles the traffic jam formation that causes car density corrugations. Besides, inhomogeneous mixing of generalized PV and a feedback process via a Rhines scale dependent mixing length were suggested as the staircases' formation mechanism in a reduced model [13,14] based on the 2D Hasegawa-Wakatani (HW) system of equations for drift-wave turbulence. Various theories and models proposed different zonal flow staircase generation mechanisms [15][16][17][18]. According to our knowledge, there is so far no clear answer as to how the staircase-like corrugations are generated in tokamak plasmas. The diversity of proposed zonal flow staircase generation mechanisms reflects the complexity of the problem, as well as the necessity and urgency to study the problem from either numerical or realistic experiments. Moreover, besides the generation mechanism, the saturation mechanism of the staircase instability is another important conundrum that has not been touched yet.
In this article, we demonstrate the generation of E × B staircase-like corrugations due to a nonlinear time delay with the first principle, state-of-the-art gyrokinetic numerical experiments of tokamak plasmas. The saturation mechanism of the staircase instability is also addressed. Numerical findings and experimental measurements in a KSTAR L-mode plasma [9] agree well with each other. The numerical platform we use is the global δ f particle-in-cell gyrokinetic plasma simulation program (gKPSP code), which solves the nonlinear gyrokinetic equations [19] for ions and bounce-averaged kinetics for trapped electrons [20]. This program is suitable for ion temperature gradient (ITG) and trapped electron mode (TEM) driven microturbulence and transport in tokamaks [21,22]. Simulation parameters are based on the 'cyclone' DIII-D H-mode shot (#81499) [23]. The normalized electron temperature gradient R/L Te ∈ [4,12], the ITG R/L Ti = 2.2, the ion (electron) density gradient R/L ni (R/L ne ) varies from 2.2 to 5. Here R is the major radius and L T(n) * = −(dln T(n) * /dr) −1 with * representing electron (e) and ion (i). At the middle radius, the inverse aspect ratio ε = r/R = 0.18, safety factor q = 1.4, magnetic shearŝ = q r/q = 0.78, T i = T e . Hydrogen is the main ion m i /m e = 1836 and plasma elongation κ = 2.0. This plasma is dominated by collisionless trapped electron mode (CTEM) microturbulence and associated transport [24,25].

Formation mechanism
Staircase-like corrugations in the mean E × B flow, strongly correlated mean electron density (gradient) and mean electron particle flux profiles have been observed in this TEM driven turbulence and transport, as shown in figure 1. In the figure, plots on the left column are time-averaged (denoted by . . . ) zonal flow V E×B , normalized electron density gradient − R/L ne and mean electron particle flux Γ e from (a) to (c). The time window for the time average is during the nonlinear saturation stage t ∈ [50, 300]R/V Ti . Here Long-lived, quasistationary staircase-like radial corrugations can be observed in those profiles with highly correlated radial locations of peaks and valleys. This implies the intrinsic causality among those physical quantities that we are willing to uncover in this article. Staircases propagate inward with a dominant low frequency ω ∼ 0.02V Ti /R and radial wave number k r ρ i ∼ 0.78, as shown in figure 1. The right column in the figure are contour plots of V E×B , R/L ne and Γ e in the (k r , ω) domain from (d) -(f ). One can estimate the staircase scale Δ ∼ 2π/k r ≈ 8ρ i , which is consistent, by measuring the interval between two corrugation peaks, given the minor radius a ∼ 220ρ i . Considering the spectrum of finite k r , the staircase-like represents the time evolution of staircase intensity in electron particle flux. In the same way, the staircase intensity of electron density δn e (t), zonal flow V E×B (t) and turbulence potential δφ(t) can be defined. Figure 2 shows the time evolution of δn e (t) (black-solid line), Γ e (t) (red-dash line) and turbulence potential δφ(t).
A staircase instability seed is initiated simultaneously up to around t = 30R/V Ti due to the nonlinear excitation of microturbulence driven by CTEM instability. The turbulence potential driven by pressure gradients generates the initial Both are normalized to their own maximum. disturbance in particle flux. The seed stays at a low level at the initial phase, and starts with a leading growth of Γ e , which then drives the growth of δn e and δφ with an evident time delay, as shown by two small arrows in the blue-dashed box. The seed eventually grows to staircase-like corrugations and starts to saturate when it reaches a large enough amplitude. The finite time delay between Γ e and δn e ( δφ) indicates that corrugations generated in Γ e cause corrugations to occur in δn e ( δφ) after a finite time delay. We note that the density δn e and potential δφ staircases grow nearly simultaneously. This implies that the  We note that the above process is interesting as it introduces us a novel mechanism to destabilize instabilities. Instability can grow from perturbations in transport first and then influence the background turbulence. Radial perturbations in electron particle flux Γ e causes corrugations in the density profile δn e , which drives similar corrugations in turbulence potential δφ and mean E × B flow through radial force balance as well as microturbulence. The E × B staircase plays a significant role as it shears the transport, a way to transfer energy to the staircase instability. This forms a feedback loop that generates staircase instability.
This feedback loop assumption that was made in reference [12] is then verified by simulations. In numerical experiments, we can break the feedback loop by artificially removing E × B flow, as indicated by sketch figure 3(a). After the feedback loop is broken, staircase-like corrugations cannot be generated in either electron particle flux or electron density profile, while an instability seed with small amplitude perturbation can be observed without a particular scale. Figure 3 presents the comparison of time-averaged radial profiles of electron density gradient (b) and electron particle flux (c) with (black line) and without (red line) zonal flow staircase. Therefore we conclude the feedback loop is necessary for the staircase-like corrugations formation. The reduced model based on the 2D HW equations demonstrates weak zonal flow shearing (vorticity) effects [14] but strong gradient of shearing effects by the deep-learning data method [26] in such a feedback loop.

Saturation mechanism
The E × B flow staircase not only shears the transport, but also saturates the staircase instability. Predator-prey [27] like behavior can be observed in the time evolution plots (figure 4(a)) in the saturation period (∼ t > 40R/V Ti ), as shown by red and black arrows in the figure. Limit cycle oscillation (LCO) [28,29] can be observed in figure 4(b) and (c), which are plots of Γ e versus V E×B for t < 70R/V Ti and t > 70R/V Ti , respectively. This predator-prey like LCO behavior between Γ e and V E×B clearly demonstrates that E × B flow staircase also saturates the instability.

Characteristics of the nonlinear time delay
One can understand the time delay by a common phenomenon in daily life, a traffic jam. In traffic transportation, there is a response time when an individual driver adjusts speed to the steady-state speed [5,30,31]. In plasma physics, density is particle number per volume. Transport of individual particles and their collective effects would therefore influence the density profile evolution. The above finding demonstrates the corrugations occur first in the particle flux, while an individual particle needs a finite time to respond to the flux corrugations, leading to the density staircases after a time delay σ. The nonlinear time delay observed in simulations can be regarded as the global collective effect of all particles' individual response time. Similarities between staircase formation and a traffic jam were introduced in a heuristic model [12], and our results are consistent with the concept. The particle density is similar to the car density; the particle transport is like traffic flow; corrugations correspond to traffic jams. The response time of an individual particle's transport is the key connecting to the adaptation time of an individual driver's response to the ambient traffic speed. Comparisons are summarized in table 1.
In order to quantitatively characterize the time delay, we adopt the time delay cross correlation statistics method, which has ever been used to detect the causality between time series by identifying the maximum in correlation with a finite time delay [32,33]. The cross correlation between Γ e and δn e with a time delay σ, denoted as δn e (t + σ) * Γ e (t), is computed. The time delay σ = −2.8 is measured at the maximum cross correlation. The maximum cross correlation coefficient is ∼0.91, which indicates a highly positive correlation between Γ e and δn e with a time delay σ = −2.8. This demonstrates that staircase-like corrugations in the electron particle transport cause similar corrugations in the electron density profile after a time σ = 2.8 delay in this case.
The time delay decreases with increasing mean flux, and corrugations are more arduous to form. Figure 5 shows the time delay σ becomes shorter as the mean particle flux diffusivity D increases. σ is calculated by the time delay cross correlation method. Here we define a quantity A ≡ R/L Te − R/L c Te characterizing the proximity to marginality  with R/L c Te 3 being the approximate linear critical electron temperature gradient for this case. We find the σ scales approximately σ ∝ A −0.67 [34]. It is necessary to point out that the time delay σ ∼ [1.6, 5]R/V Ti is slightly larger than the turbulence autocorrelation time ∼ [1,3]R/V Ti . Figure 6 presents the formation of staircase-like corrugations in three cases from low to high particle flux (top to bottom). It can be seen from the figure that profile corrugations disappear by increasing the particle flux. This demonstrates that there should be a critical time delay σ c . When the time delay σ < σ c by increasing the mean particle flux, corrugations cannot be generated. The numerical experiments of various parameters yield self-consistent characteristics of the staircase-like pattern formation.

Supporting evidence from experiments
On the other hand, a KSTAR experiment of L-mode plasmas yields consistent results with the numerical findings. This plasma is MHD instability free and characterized by the co-existence of avalanche-like events and electron temperature profile corrugations (or E × B staircases). More details about the experiments can be found in reference [9]. We revisit the experimental measurements, and the time evolution of electron temperature gradient R/L Te near the avalanche initiation location and vertical profiles of electron temperature variation δT e (Z) = T e − T e across the plasma center (Z = 0) are shown in figure 7. Here, T e means the long time average of T e for 40 ms. The figure presents several small and large avalanche events characterized by a sudden drop of the gradient for the selected time range t = 3.9 s to t = 3.93 s (some of them are indicated by arrows in figure 7). After an avalanche events occur, the electron temperature gradient starts recovering and increasing slowly. Vertical profiles of electron temperature variation δT e (Z) are plotted for several particular time points t 1 = 3.90666s, t 2 = 3.90879s, t 3 = 3.90982s, t 4 = 3.91225s, t 5 = 3.9176s, t 6 = 3.9194s.
Understanding the mechanism of formation and destruction of T e corrugation (or the E × B staircases) is important, since it can limit the avalanche events and influence transport barrier formation. T e corrugations observed around two time points (t = t 1 and t 4 ) are shown for examples. When the gradient reaches a high enough value exceeding some threshold and the corrugation disappears (eg. t = t 2 and t 5 ), a large scale heat avalanche-like event can occur [9]. Destruction of the T e corrugation with large T e gradient is understandable from the aforementioned simulation results. The electron temperature gradient enhances the background turbulent heat flux and therefore reduces the time delay σ. When the σ is shorter than the critical value σ < σ c , the staircases will be destroyed. The temperature corrugation can reform when the gradient is reduced sufficiently and the heat flux to be jammed is provided by avalanche like events.

Conclusions and discussions
In summary, numerical experiments show that there exists a finite time delay between the staircase-like corrugations formation in particle transport and density profile (turbulence potential). This demonstrates a novel picture that instability starts from perturbations in transport and then influences the background turbulence. E × B flow staircase shears the transport, which is necessary to transfer energy in the feedback loop of the staircase generation. On the other hand, the E × B flow staircase also plays a significant role as the main saturation mechanism for the staircase formation. Findings from numerical experiments are consistent with realistic experiments in KSTAR tokamak L-mode plasmas.
Finally, based on the numerical and experimental results, we can propose a dynamic process loop among profile gradient, zonal flow staircase and some time delay, as shown in figure 8. We note that this dynamic process loop can be a fresh heuristic mechanism to understand the edge transport barrier formation related phenomena, such as periodic L-H-L transitions [35], intermediate state (I phase) [36] during L-H transition [37] and edge localized mode (ELM) cycles [38] etc.