{'fraction_non_alphanumeric': 0.0313204146011717, 'fraction_numerical': 0.013970256872465074, 'mean_word_length': 4.830213464696223, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Saturation of Kinetic Ballooning Modes by Zonal fields in Toroidal Plasmas

Kinetic ballooning modes (KBM) are widely believed to play a critical role in disruptive dynamics as well as turbulent transport in tokamaks. While the nonlinear evolution of ballooning modes has been proposed as a mechanism for detonation in tokamak plasmas, the role of kinetic effects in such nonlinear dynamics remains largely unexplored. In this work global gyrokinetic simulation results of KBM nonlinear behavior are presented. Instead of the finite-time singularity predicted by ideal MHD theory, the kinetic instability is shown to develop into an intermediate nonlinear regime of exponential growth, followed by a nonlinear saturation regulated by spontaneously generated zonal fields. In the intermediate nonlinear regime, rapid growth of localized current sheet is observed.

Ballooning instability (or its astrophysical counterpart, the Parker instability) in a magnetized plasma is characterized by long wavelength parallel and short wavelength perpendicular to the background magnetic field. The nonlinear evolution of the instability has been a subject of great interest for a diverse range of eruptive phenomena in space plasmas such as substorms in the Earth's magnetotail [1,2] and edge-localized modes (ELMs) [3] and disruption precursor modes [4] in toroidal fusion plasmas. Theoretical studies of nonlinear ideal magnetohydrodynamic (MHD) ballooning modes [5,6] predict explosive nonlinear growth near the linear instability threshold. Finger-like structures develop, forming a front with steep pressure gradient which can further destabilize the mode nonlinearly, and result in a finite-time singularity ("detonation") [5]. However, attempts at simulating such an instability using the full MHD equations have not succeeded in realizing the tendency to form a finite-time singularity. While finger-like structures and the steepening of pressure gradients are indeed observed [7] as predicted by the analytical theory, the mode is seen to grow in the nonlinear regime exponentially with approximately its linear growth rate, maintaining its filamentary structure. A new asymptotic regime, called the "intermediate" nonlinear regime of exponential growth has been formulated analytically to account for these simulations [8].
Exponential nonlinear growth, even though it lacks the impulsive signature of "detonation", is of great interest. During this phase, the width of the finger-like structures becomes sufficiently small that the validity of the MHD model is questionable. In weakly collisional plasmas, kinetic effects will intervene. This leads to considerations of the kinetic ballooning mode (KBM) which is widely recognized to play an important role in fusion [8,9] as well as space plasmas [10]. However, the nonlinear dynamics of the KBM, and in particular the role of kinetic effects [11] is not well understood. Numerical simulations on the nonlinear saturation of KBMs have been inconclusive [12][13][14]. Various simulation results suggest that KBMs can only saturate nonlinearly with pressure profile relaxation [12] or increased flow shear [13] above a critical value of beta.
In this Letter, we demonstrate from a gyrokinetic particle-in-cell (PIC) global simulation that following a linear regime of exponential growth, the KBM evolves into an intermediate regime which exhibits features that are qualitatively different from its ideal MHD counterpart [7]. First, the growth is slightly faster than linear exponential growth, which indicates that the perfect cancellation that occurs in the intermediate ideal MHD dynamics between nonlinear destabilization due to enhanced pressure gradients and stabilization due to field-line bending [8] does not occur in the kinetic intermediate regime. Second, and more important than the first, the kinetic electromagnetic dynamics lead to the spontaneous generation of zonal flow (flux-surfaceaveraged electrostatic potential δφ ) and zonal current (flux-surface-averaged vector potential δA ) through three-wave coupling processes, in contrast to the generation of zonal flow through modulational instability in electrostatic ion-temperature-gradient (ITG) turbulence [15]. When the zonal flow shear exceeds the linear growth rate, zonal flow shearing suppresses the nonlinear instability which in turn, self-regulates the zonal fields (the zonal flow and the zonal current), leading to a saturated nonlinear regime. In the kinetic intermediate regime, thin current sheets develop near the mode rational surfaces, which can eventually exhibit tearing instability in the presence of resistivity, but the resistive taring mode growth rate appears to be too slow to have a strong effect on the nonlinear dynamics of the instability, which is dominated by the regulation of the mode by zonal fields.
Gyrokinetic simulation of KBM.-In the simulations using the gyrokinetic toroidal code (GTC) [16,17], ions are treated by the gyrokinetic Vlasov equation, while electrons are described using the nonlinear fluid equations: the electron perturbed density δn e is calculated by time-advancing the continuity equation [18], and the electron parallel flow δu e is calculated by inverting the parallel Ampere's Law [19]. The gyrokinetic Poisson's equation is solved to obtain the perturbed electrostatic potential δφ. The parallel vector potential δA = δA adi +δA na is solved for the adiabatic and non-adiabatic parts. Integrating the electron drift kinetic equation to the momentum order and considering the parallel Ampere's law, we arXiv:1805.00402v1 [physics.plasm-ph] 1 May 2018 can derive the linear Ohm's law for adiabatic δA adi and the nonlinear Ohm's law for non-adiabatic δA na as: where B 0 is the equilibrium magnetic field, δB ⊥ = δB adi ⊥ + δB na ⊥ is the perturbed perpendicular magnetic field, which consists of the adiabatic and non-adiabatic parts as δB adi and δB is the compressional magnetic perturbation, which is solved using the gyrokinetic perpendicular force balance equation [21]. Here ∂ψ0 −δφ is the inductive potential. δψ adi is the adiabatic component of the perturbed poloidal flux, which is defined as ∇δψ adi × ∇α = δB adi ⊥ , where α = q(ψ 0 )θ−ζ is the fieldline label with the Boozer poloidal angle θ and toroidal angle ζ, and the safety factor q(ψ 0 ) is a function of the equilibrium poloidal flux ψ 0 . Also, T e is the electron equilibrium temperature, n 0 is the plasma equilibrium density, and P e0 = n 0 T e is the electron equilibrium pressure. The first term on the right-hand-side of Eq. (2) represents the nonlinear ponderomotive force, which is shown in the following sections to play an important role in KBM nonlinear dynamics. The second and third terms are the nonlinear drive from finite δB , and are small compared with the nonlinear ponderomotive drive due to the smallness of β. A complete form of the generalized Ohm's law is presented in [20]. In addition to the zonal current, the nonlinear ponderomotive force in Eq. (2) can produce localized current sheets on the mode rational surfaces with tearing parity. Since the typical current sheet width realized in our simulations is larger than the electron skin depth, we drop the electron inertia terms in the generalized Ohm's law. In the simulations, Cyclone Base Case parameters are used for the background plasmas: the major radius is R 0 = 83.5cm, inverse aspect ratio is a/R 0 = 0.357. At r = 0.5a, the plasma parameters are B 0 = 2.01T , T e = 2223eV , R 0 /L T = 6.9, R 0 /L n = 2.2, q = 1.4, β e = 2%. The first order s − α model with circular cross-section [18] is used for the equilibrium magnetic field, with θ = θ 0 − sin θ 0 + O( 2 ), where θ 0 is the geometric poloidal angle, and = r/R 0 is the normalized radial coordinate. We simulate n=10 toroidal mode KBM (keeping all the poloidal harmonics m), and its nonlinear interaction with the zonal mode (m=0, n=0). In the linear simulations, the mode exhibits ballooning mode structures, with a real frequency ω lin r = 0.67c s /a and maximum growth rate γ lin = 0.62c s /a. Convergence studies show that the physical results in the linear and nonlinear simulations are not sensitive to grid size, time step size or number of particles per cell.
Saturation by zonal fields.-A time history for the nonlinear KBM simulation is shown in Figure 1. The perturbed electrostatic potential, parallel vector potential, and parallel magnetic field are normalized, consistent with ideal MHD scaling, as eδφ/T e , cδA /v A B 0 R 0 and δB /B 0 respectively, where v A = B 0 / √ 4πn 0 m i . Normalization will be omitted in the following text for simplicity. The perturbed electrostatic potential δφ 10,14 , the parallel vector potential δA 10,14 , and the parallel magnetic field δB 10,14 of the dominant (10,14) mode are measured at the mode rational surface with q = 1.4 at the center of the simulation domain, and the zonal flow δφ and the zonal current δA amplitude are averaged over the simulation domain. The mode starts as a small perturbation in the fields, and quickly evolves into the linear regime after a brief transient stage. Before t ∼ 11a/c s , δφ 10,14 grew more than two orders of magnitudes at the linear growth rate γ lin . Here δA 10,14 remains around three orders of magnitude lower than δφ 10,14 , since the linear adiabatic component δA adi 10,14 is zero at the resonant surface, as constrained by Eq. (1). Also, δB is more than one order of magnitude smaller than δφ due to the small plasma β. At t ∼ 11a/c s , δA 10,14 starts to grow faster than exponential, indicating that the mode evolves into an intermediate regime, where the nonlinear ponderomotive effects become important. From t ∼ 11a/c s to t ∼ 15a/c s , δφ 10,14 and δB 10,14 grow slightly faster than exponential, with an effective growth rate γ int = 1.1γ lin between t = 11a/c s and t = 14.5a/c s , and the field quantities retain their linear poloidal mode structure. These features of KBM in this intermediate regime are similar to those in the intermediate regime found in compressible MHD simulations [7]. The growth of dominant field quantities at a rate faster than predicted by linear theory indicates that exact cancellation of nonlinearities in the ideal MHD intermediate regime does not occur in this kinetic intermediate regime [8]. In the linear regime and the intermediate regime, the zonal flow and the zonal current both grow exponentially at a growth rate γ zonal ∼ 2γ lin . This shows that the zonal fields in KBM are driven more efficiently by the threewave parametric instability, in contrast to the zonal flow excitation by modulational instabilities in electrostatic ITG, where the zonal flow grows exponential on exponential and the instantaneous growth rate depends on mode amplitude [15].
At t ∼ 15a/c s , the dominant mode and the zonal fields saturate nonlinearly. As shown by the diamond dotted blue line in Figure 1, the steady state zonal flow amplitude is around 5 times larger than the dominant δφ 10,14 component. At the same time, the ion energy transport reaches steady state at the gyro-Bohm level with χ i ∼ χ GB , as shown by the black solid line in Figure 2, where δf is defined as the volume averaged ion energy flux normalized by local temperature gradient, where v r is the radial drift velocity including the E × B drift and the magnetic flutter drift [20]. In the simulations where the zonal flow and the zonal current are both artificially suppressed, the nonlinear ion heat conductivity becomes one order of magnitude larger, as shown by the diamond red line in Figure 2. Note that δφ 10,14 also saturates at a magnitude around 3 times higher than that in the case with both zonal flow and zonal current, as shown by the diamond dotted blue line in Figure 2. In two other simulations where only the zonal current or the zonal flow is artificially suppressed, δφ and χ i saturation levels also see a significant increase, indicating that both zonal flow and zonal current regulate ion energy transport and mode saturation in KBM. A comparison of δφ nonlinear poloidal structure between simulations with and without the zonal fields is shown in Figure 3. In the simulation with selfconsistently generated zonal flow and zonal current, the zonal fields break up the radially elongated eigenmode structure into micro and mesoscale structures, reducing radial transport. In the simulation with the zonal fields artificially suppressed, although the non-zonal nonlinear E × B term also shears the mode structure, some macroscale radial filaments of streamers survive. These results show that the KBM saturation is governed by the zonal fields, including both the zonal flow and the zonal current. In two additional simulations where β e = 1.74% and β e = 1.55% (near the KBM instability threshold), we observe similar nonlinear saturation features.
Onset of nonlinear rapid growth of localized current sheet.-As shown by the diamond solid red line in Figure 1, δA 14,10 at the mode rational surface first grows faster than exponential, and then grows more than one order of magnitude exponentially with a nonlinear growth rate γ nl ∼ 3γ lin during the intermediate nonlinear regime. The nonlinear growth rate can be explained by the coupling between the zonal current and non-zonal inductive potential through the first term in Eq. (2). The poloidal δA structure evolves from the linear eigenmode structure at t = 11c s /a, as shown in Figure 4(a), to mesoscale structures at t = 17c s /a, as shown in Figure 4(b). The mode structure becomes very thin in the radial direction, and remains the same scale as the linear eigenmode structure in the poloidal direction. This corresponds to the rapid growth of current sheets localized at the rational surfaces, excited by the nonlinear ponderomotive force terms in the generalized Ohm's law. In the simulation where the nonlinear ponderomotive force terms are not included (δA na = 0), although zonal flows still break the linear mode into mesoscale structures nearly isotropic in radial and poloidal directions, as shown in Figure 4(c), the radial correlation length of the turbulence eddies is much longer than that in the case with the self-consistent ponderomotive force.
The development of localized current sheet in the intermediate and nonlinear regime in KBM is analogous to the nonlinear process in the ideal MHD theory. However in this scenario where the kinetic effects become important during the intermediate regime, the mode saturates at the spatial scale comparable to the ion gyroradius with a transport level controlled by the zonal fields. In contrast, the mode structure in the MHD theory tends to become singular until the pressure profile flattens by transport. The radial profiles of (n,m) harmonic of δA at t = 11c s /a and t = 17c s /a are shown in Figure 4(d) and Figure 4(e). The linear mode structure has exact odd parity at the rational surfaces, and the nonlinear mode structure contains even parity component at the rational surfaces driven by the nonlinear electromagnetic ponderomotive force. For comparison, Figure 4(f) shows the (n,m) harmonic of δA after saturation in the simulation with δA na = 0. In this case, although the mode structures deviate from the linear mode structures, each (n,m) harmonic is still zero at the q = m/n surface. Finally, the mode exhibits a forward cascade in the radial mode number k r as it evolves to the nonlinear regime, corresponding to the nonlinear shearing and breaking of mode structures as shown in the poloidal contours. Because of the formation of thin current layer near rational surfaces, we conducted simulations with finite resistivity in the generalized Ohm's law to test the role of resistive tearing physics in the saturation of KBM [22]. With resistivity 100 times the Spitzer resistivity, no significant tear-ing instability is observed within the time scale of KBM nonlinear saturation. In this case, KBM linear growth rate and real frequency are increased significantly by the resistive drive, and the zonal fields still saturate the mode with radially smoother nonlinear mode structure.
Conclusions and future work.-In summary, we have presented global gyrokinetic particle-in-cell simulation results of KBM nonlinear behavior. Instead of finitetime singularity, the instability develops into an intermediate nonlinear regime of exponential growth, followed by nonlinear saturation regulated by spontaneously generated zonal fields. Both zonal flow and zonal current are important for the saturation level. In the intermediate nonlinear regime, onset of rapid nonlinear growth of localized current sheet is observed. Although we did not observe significant growth of resistive tearing mode on the time scale of KBM nonlinear evolution, the current sheet near the rational surfaces might induce tearing instabilities on a longer time scale, which can affect the transport level, and provide seed islands for global tearing modes such as the neoclassical tearing mode. Therefore including the collisionless tearing physics by keeping the electron inertia terms in the simulation model [20] is an important next step. On the other hand, KBM turbulence can induce large transport at plasma edge in H-mode plasmas [23][24][25]. According to the EPED model [26], KBM and peeling ballooning mode (which has both ballooning drive and bootstrap current drive) together determines the pedestal height and width. In future work, we plan to extend the KBM simulations presented in this Letter to experimentally measured equilibrium profiles with multiple modes to determine the prediction of KBM turbulence transport, its coupling with tearing modes and its role in shaping the edge.