A fast edge turbulence suppression event has been simulated in the electrostatic version of the gyrokinetic particle-in-cell code XGC1 in a realistic diverted tokamak edge geometry under neutral particle recycling. The results show that the sequence of turbulent Reynolds stress followed by neoclassical ion orbit-loss driven together conspire to form the sustaining radial electric field shear and to quench turbulent transport just inside the last closed magnetic flux surface. The main suppression action is located in a thin radial layer around $\psi N\u22430.96\u20130.98$, where *ψ _{N}* is the normalized poloidal flux, with the time scale $\u223c0.1$ ms.

## I. INTRODUCTION

Transport barrier formation and its relation to the flow of the fluid medium are of fundamental interest in various natural and laboratory observations, such as geophysical and atmospheric fluid systems.^{1,2} In a magnetic fusion device, this physics has a critical implication to achieving an economical energy production since the bifurcated plasma state, called the high confinement mode (H-mode),^{3} is often envisioned as the operating mode of choice for fusion reactors,^{4} and will be relied in ITER in achieving its goal of ten-fold energy gain.^{5} However, despite over 30 years of routine H-mode operation in all the major tokamaks, there has been no fundamental understanding at the kinetic level on how the H-mode bifurcation occurs.

Experimentally, a radial transport bifurcation into the H-mode in both plasma density and thermal channels occurs in a thin edge layer of the tokamak plasma just inside the magnetic separatrix surface when the plasma heating power exceeds a critical value.^{3} As a result, a plasma density and temperature pedestal is formed on the time scale of a few ms with a steep gradient in the thin edge layer. As this pedestal forms, the core plasma pressure inside the edge layer position increases, resulting in a transition of plasma operation to a high confinement H-mode from a low confinement L-mode.^{3} The bifurcation event is accompanied by a sharp increase in the sheared E × B flow and a significant drop in the turbulence amplitude within the thin transport barrier layer on a time-scale that is often shorter than 0.1 ms if the heating power is strong (strongly driven). The edge heating needed to initiate this H-mode regime is minimized when the ion magnetic-drift direction, or the $B\u2192\xd7\u2207B$-direction, is toward the magnetic X-point when the plasma is operated with a single poloidal divertor.^{6}

There have been many attempts to apply simple theoretical models to explain how an H-mode transition could occur. A popular “predator-prey” model^{7} implies that increasing the heat flux to the edge of the plasma, thus raising the edge gradients, results in stronger turbulence (prey). The increased turbulence could then amplify the sheared poloidal flow (predator) nonlinearly through the turbulent Reynolds stress. When the flow drive is larger than the flow damping, the sheared poloidal flow could grow, nonlinearly extracting even more kinetic energy from the turbulence. As a result, the turbulence and the associated turbulent transport could collapse. This suppressed turbulence state is then conjectured to be maintained through the steep-pressure driven sheared E × B flow driven by the simultaneous build-up of the H-mode pedestal.

Extended predator-prey models predict both an oscillatory limit-cycle (LCO) type predator-prey transition^{8} and a sharp transition^{9–11} triggered by a single burst of axisymmetric sheared turbulence-driven E × B flow (known as zonal flow). Turbulent fluid simulations have shown evidence for some of this phenomenology.^{12–14} Laboratory experiments have indeed reported both LCO type transition^{15–20} when operating close to the H-mode power threshold, and a sharp bifurcation^{21–24} within 0.1 ms (Ref. 24) when the power threshold is exceeded. In the fast transition, some detailed experiments report that the turbulent stress-driven shear flow first leads to a collapse of the turbulence, which is then followed by the development of the edge pedestal in a rather longer time scale, claiming that the turbulence suppression is not maintained by the simultaneous buildup of the steep pedestal and the associated E × B shearing.^{25,26} Some experiments (i.e., Ref. 27) report a different evidence that the experimentally observed Reynolds work is too weak to explain the L-H bifurcation and, thus, the E × B shearing from the neoclassical orbit loss physics^{28,29} is solely responsible for the bifurcation. Some LCO type transition experiments that showed little Reynolds work suggest that the E × B shearing from the buildup of the steep pressure pedestal is responsible for the turbulence suppression.^{30}

This body of evidence suggests that the H-mode transition could indeed be related to the sheared E × B-flow, either turbulence-, orbit-loss- or $\u2207p$-driven. However, the existing models are based upon simplified ad-hoc equations and the turbulence simulations assume specific instability mechanisms, ignore possible important kinetic effects, and are not carried out in a realistic geometry. This paper presents a more detailed report on the first study of edge transport barrier formation dynamics in the Physical Review Letter article,^{31} using a first-principles based electrostatic gyrokinetic simulation implemented in XGC1 in realistic edge geometry.^{32,33} In the gyrokinetic equations, the fast gyro-motion is analytically treated, thus removing the gyrophase angle variable, while preserving the most basic plasma physics element at first principles level; i.e., the motions of individual particles and their parallel Landau resonance with waves. Moreover, the XGC1 simulations evolve the total distribution functions $f(x,v,t)$ for the gyrokinetic ions and drift-kinetic electrons, hence the background macro-scale kinetic neoclassical physics is self-consistently included together with the micro-scale nonlinear turbulence physics and no *a priori* instability-drive assumption is made, except the low-beta electrostatic-limit assumption. Even though the boundary plasmas in most of the tokamak experiments are at electron *β _{e}* (=electron kinetic energy/equilibrium magnetic field energy) close to or below the electromagnetic criticality (= $me/mi\u22430.03%$ for deuterons), they may be influenced by electromagnetic effect through some modification of the parallel electric field fluctuations. The study of the electromagnetic effect on the boundary plasma will be reported elsewhere in the future. In the present simulation of the Alcator C-Mod

^{34}L-mode plasma discharge #1140613017, the plasma

*β*in the relevant edge region is only $\u2243$0.01%, which is well below the electromagnetic criticality. Thus, we assume that the electromagnetic effect on the edge turbulence suppression physics is much weaker than the electrostatic effect. This assumption needs to be verified in the future. Experimental observations of the “blobby” boundary turbulence have shown evidence for dominance of electrostatic fluctuations.

_{e}^{35}In order to handle the orbit-loss hole and non-Maxwellian physics properly in the velocity space, a conserving and fully nonlinear Fokker-Planck collision operator is used.

^{36}Lost plasma particles are recycled as Monte Carlo neutral atoms in the divertor chamber, with charge exchange and ionization interactions with plasma.

To enable the difficult boundary plasma simulation in a gyrokinetic code, a new multi-scale, total-*δf* particle-simulation scheme called a “hybrid-Lagrangian” scheme has been developed and used in XGC.^{32} The scheme utilizes a velocity-space grid in handling a non-Maxwellian plasma, in performing a fully nonlinear Fokker-Planck collision operation, in reducing particle noise, and other boundary-specific issues to be described in this paper. In Ref. 32, the scheme has been used to verify its effectiveness in the simulation of the ion temperature gradient (ITG) turbulence with gyrokinetic ions and adiabatic electrons. In the boundary plasma, a proper kinetic electron scheme is essential since first-principles-based dynamics of kinetic electrons is critical in producing the first-principles-based 3-dimensional mean electric field vector, the sheared E × B flow, and the proper micro-turbulence modes in the boundary plasma consistently with each other. Moreover, the plasma ions and electrons are lost to the material wall, which are then recycled back into the plasma with atomic charge exchange and ionization interactions. Under the wall-loss conditions, the phase-space volume conservation requirement is a non-trivial issue. At the same time, a gyrokinetic-consistent sheath boundary condition should be used to keep the plasma quasi-neutral. It will be shown in the present report how the hybrid-Lagrangian scheme can be utilized to resolve these difficult issues.

Kinetic electrons not only add physics difficulty when in contact with the material wall, but also high computational cost and algorithmic complexity in the computer hardware accelerators such as Graphics Processing Units (GPU). The electron parallel speed that is faster than that of the ions by a factor $ma/me$ ( $\u224860$ for deuterons), with *m _{a}* being the mass of species

*a*, makes the simulation time-step smaller by the same factor. The large amplitude turbulent fluctuations in the boundary plasma could trigger a more sensitive Courant instability issue.

^{37}Due to the fast parallel motion of the tail electrons in the direction parallel to the magnetic field, domain decomposition is also difficult. In order to address these difficulties, special schemes are needed.

Therefore, we take a progressive approach in attacking the difficult boundary plasma physics and solve the electrostatic problem first, with the electromagnetic turbulence corrections to be investigated afterwards. The electrostatic problem can be simpler in some aspects, but it introduces an additional complication due to the so called “*ω _{H}* mode” instability

^{38}from the catastrophic interaction of the passing kinetic electrons with the electrostatic-version of the shear-Alfven wave. The edge physics model and algorithm must be able to cope with this instability to obtain robust solutions. We study the low-recycling regime first, leaving the study of the high-recycling regime to a later time with more accurate atomic physics interactions that are adequate for lower divertor temperatures, such as molecular physics and volume recombination.

This paper is organized as follows: in Sec. II, we briefly summarize the electrostatic version XGC1 and the gyrokinetic equations,^{39} and the basic features of the hybrid-Lagrangian scheme that utilizes a 2-dimensional (2D) velocity space grid at each particle-in-cell grid node.^{32} In Sec. III, we describe the unusual physics issues that are specific to the boundary plasma, followed by the new physics models and the enabling algorithms adopted in XGC1 by taking advantage of the hybrid-Lagrangian scheme. We then verify the kinetic electron dynamics in the hybrid-Lagrangian scheme by cross-comparing the ITG-TEM modes transition against the published results, and demonstrate the scheme's capability in generating the nonlinear turbulence across the magnetic separatrix. Here, “ITG modes” stands for the “ion temperature gradient modes” and “TEM” stands for the “trapped electron modes.” In Sec. IV, we describe the simulation setup and present the dynamics of the edge turbulence suppression in more detail than the abbreviated report in the Physical Review Letter article.^{31} The conclusions and discussions are summarized in Sec. V. XGCa is the axisymmetric version of XGC1 for gyrokinetic neoclassical physics, uses the same equations and the hybrid-Lagrangian scheme, and will not be specifically discussed in the present report.

## II. A BRIEF INTRODUCTION TO THE XGC1 CODE AND THE BASICS OF THE HYBRID-LAGRANGIAN SCHEME

The material presented in this section is not new and published elsewhere as cited, but is summarized here for completeness sake in providing the readers a self-contained background information before progressing to the new and more advanced algorithms for the boundary plasma in Sec. III.

XGC1 is a multiscale 5-dimensional (5D) gyrokinetic code that solves the turbulence-neoclassical-neutral transport physics together, and is specialized for tokamak edge plasma based on the particle-in-cell technique combined with a velocity-space grid technique (hybrid-Lagrangian scheme) in the most recent version described here. The simulation domain normally includes the whole plasma volume—from the material wall, across the magnetic separatrix surface with the X-point, and to the magnetic axis—in order to avoid artificial Dirichlet or Neumann boundary conditions at the core side. Unstructured triangular meshes that follow approximately the equilibrium magnetic field lines are used to describe the complicated wall geometry and the X-point.

The only solver boundary condition is the grounded zero electrostatic potential at the material wall, which is approximated as axisymmetric. Once a magnetic field-line intersects a material limiter other than in the divertor chamber, it is assumed that the strong electrical conductivity holds the electrostatic potential to be zero along the magnetic field lines. This constraint on the solver boundary condition does not apply to the plasma or neutral particle motions: they see the real material wall geometry. The material wall geometry, the magnetic equilibrium profiles, and the initial plasma profiles are imported from experimental data using the magnetohydrodynamics equilibrium fitting code, EFIT^{40} outputs, or from an analytic Grad-Shafranov equilibrium. Wall loss of the plasma particles induces the birth of Monte Carlo neutral atoms at the Frank-Condon energy at a specified recycling rate.

### A. Gyrokinetic equations

The electrostatic version XGC1 solves the 5D gyrokinetic Boltzmann equation^{41,42} in the following form:

Here, *f* is the distribution function of gyrokinetic particles, $X$ is the gyro-center position in the configuration space, *S*(*f*) represents the dissipative terms (such as Coulomb collisions, heating/cooling sources, and atomic interactions with neutral particles), $v||$ is the velocity of the gyro-center parallel to the local magnetic field $B,\u2009b\u0302=B/B$, $\mu =mv\u22a52/2B$ is the magnetic moment, $v\u22a5$ is the perpendicular velocity to the local magnetic field vector, $E\xaf$ is the gyro-averaged electric field, *m* is the mass, and *q* is the particle charge.

In the present study, the electrostatic potential is determined by the lowest-order quasi-neutrality equation that is also called the “gyrokinetic Poisson equation”^{42}

where $n\xafi$ is the ion gyro-center density (not the real ion density) averaged over the gyro-orbit, *n _{e}* is the electron density, $\u2207\u22a5$ is the gradient operator perpendicular to the magnetic field vector $B$, and

*ρ*is thermal ion gyroradius. For the sake of the discussion in the present paper, we use hydrogenic atoms and replace

_{i}*q*by the unit elementary charge “e” from this point on. XGC1 actually has multiple impurity species capability. Electrons are drift kinetic and the gyro-averaged electric field $E\xaf$ in Eq. (1) is replaced by $E$ in the electron equation of motion, unless the turbulence contains the electron scale modes, such as the electron temperature gradient driven modes (ETG modes).

^{43,44}

### B. The basics of the hybrid-Lagrangian total-*δf* scheme

Again, the summary presented in this subsection is not new, but is for completeness sake. This is a necessary basic material in understanding the more advanced scheme presented in Sec. III.

In solving the gyrokinetic Boltzmann equation, the most recent version XGC1 uses the hybrid-Lagrangian total-*δf* scheme,^{32} which utilizes a 2D (2-dimensional) structured rectangular grid in the velocity space in addition to 3D unstructured triangular particle-in-cell grids in the configuration space. If we define *D*/*Dt* as the left hand side of Eq. (1), the equation can be represented by

Note that *D*/*Dt* operation itself must conserve the phase space volume^{42} when the dissipation source *S*(*f*) is absent.

In the present hybrid-Lagrangian total-*δf* scheme, the distribution function *f* is decomposed into $f=f0+\delta f$, where $\delta f=fp$ is described by marker particles and $f0=fa+fg$, with *f _{a}* being an analytic function and

*f*residing in the 5D grid-space. Equation (3) can then be equivalently written in the form

_{g}If *f*_{0} is set to zero, XGC1 becomes a pure “full-f” particle-in-cell code. If *f _{a}* is set to be a time-invariant, flux-function Maxwellian

*f*with

_{M}*f*= 0 and the operator $Df0/Dt$ on the right hand side is simplified to include only the fluctuation-driven term (by neglecting the magnetic $\u2207$B and curvature-drift driven terms), XGC1 becomes a so-called “reduced

_{g}*δf*” code with $\delta f=fp$. In the reduced

*δf*formalism, a steep radial gradient in

*f*is problematic since it makes the magnetic drift term non-negligible compared to the fluctuation driven term on the right hand side due to the neoclassical orbit excursion. At the same time, if a linearized Fokker-Planck collision operator were used, the condition $\delta f\u226afM$ would be required. These simplifications cannot be enforced in the edge plasma where the magnitude of

_{M}*δf*can be $O(fM)$, and the radial banana drift width by the magnetic drift is of the same order as the radial-gradient scale length. Thus, XGC1 assumes neither that the distribution function $f=fa+fg+fp$ is near-Maxwellian nor that the magnetic drift in the $Df0/Dt$ driving-force is negligible. It has been shown that this total-

*δf*scheme yields the same self-organized quasi-equilibrium state as the full-f scheme does,

^{39}even though the transient behavior may be different due to different noise level. As $\delta f=fp$ grows, a small fraction of it can be transferred to

*f*at every time step time to mitigate the usual “growing-weight problem,”

_{g}^{45}as demonstrated in Ref. 32. The Maxwellian part of

*f*can also be transferred to

_{g}*f*while preserving

_{M}*f*

_{0}. The details can be found in Ref. 32.

In order to describe the effect of the *S*(*f*) term and the steep density gradient profile in the edge plasma on the particle distribution function *f _{p}*, we use the two-weight scheme.

^{32,46,47}For simplicity, we choose the initial

*f*to be

*f*, i.e., $fg(t=0)$ = $fp(t=0)=0$. The full-f weight

_{a}*w*

_{0}is defined as $w0g\u2261f0(zt=0,t=0)$, where $zt=0$ represents initial positions of the marker particles in the 5D phase space and

*g*is the spatially uniform distribution function of the marker particles. The initial

*w*

_{0}defines the density gradient.

*w*

_{0}is invariant in the dissipationless time-advance of the marker particles, i.e

*.,*the operation on the left-hand-side of Eq. (4).

*f*is then the time-advancing marker particle distribution function with the time-advancing

_{p}*δf*weight

*w*

_{1}, $fp=w1w0g=w1f0(zt=0,t=0)$. $w1=0$ at $t=0.$

Particle weight evolution from the driver $Df0/Dt$ is evaluated using the direct weight evolution method in the difference form.^{32} The usual differential form of weight evolution can be written as

Note that *g* is constant along a collisionless and sourceless Vlasov particle trajectory since the phase space volume is conserved. This causes a problem when the marker particles are lost at the material surface. A special technique will be introduced in Sec. III D to remedy this problem. Time integration of Eq. (5) gives a simple difference form of weight change

where Δ represents the finite difference along the particle trajectory. The Euler method can be applied to obtain the time integration of *S*(*f*). A key advantage of using the difference form is that it can avoid a numerical derivative of *f*_{0}. This is important especially when we use a noisy interpolated quantity *f _{g}* on the 5D phase space grid. Again, more details can be found in Ref. 32.

In this scheme, the entire *f* can be presented on the 5D continuum grid by mapping *f _{a}* and

*f*to the velocity grid and adding to

_{p}*f*.

_{g}*S*(

*f*) of Eq. (4) is then evaluated on the 5D continuum grid and its effect on

*f*can be transferred back to particle weights using the inverse mapping. If a marker particle corresponding to the 5D continuum grid does not exist, the

_{p}*f*-information can be left on the velocity grid. For verification of the nonlinear Coulomb operation using this scheme, we refer the readers to Refs. 36 and 48. A description on using this scheme for the atomic interaction with neutral particles is presented in Sec. III G.

## III. UNUSUAL PHYSICS CONDITIONS, PHYSICS MODELS, AND THE EXTENDED HYBRID LAGRANGIAN SCHEME IN THE BOUNDARY PLASMA

Boundary plasma physics has extra features that are absent from conventional core plasma physics, hence it requires special treatment. First of all, the plasma particles that hit the material wall are lost from the simulation domain. The requirement for the phase-space volume conservation of the Vlasov system (time-invariance of the function *g*) can be broken if a special treatment is not employed. Second, the lost plasma particles come back into the plasma in the form of neutral atoms after the molecular dissociation into atoms (we presently neglect in XGC1 the molecular dissociation procedures in front of the divertor plates in the present study of low-recycling divertor plasma). The neutral atoms then interact with the plasma particles through ionization and charge-exchange processes. Some of the plasma and neutral particles even escape the system via vacuum pumping and material binding. Third, the open system also drives the plasma to a non-Maxwellian steady state, nullifying the convenient orderings used for the core plasma that originated from the near-thermal-equilibrium assumptions. One of the consequences is that the Coulomb collisions need to be described by the fully nonlinear Fokker-Planck operator. Development and utilization in the XGC family codes of a fully nonlinear Fokker-Planck collision operator was already reported in Refs. 36 and 48, and will not be a subject of discussion here.

### A. Electron subcycling

The wave frequency *ω* of interest in XGC1 is much lower than $k\u2225ve,th$, where $k\u2225$ is the wave vector component in the magnetic field direction and $ve,th$ is the electron thermal speed. At such frequencies, there is a significant adiabatic electron response to the electrostatic potential perturbation on a flux surface. In an edge plasma, both the mean and fluctuating electrostatic potential variation on a flux surface can be significant. To reduce the particle noise in the adiabatic response to the electrostatic potential perturbation, the analytic part of the electron distribution function $fae$ is chosen to be Maxwell-Boltzmann with the perturbed potential containing both the mean and fluctuating parts as they develop during the simulation,^{49}

where *n*_{0} is the initial flux-function equilibrium electron density that is equal to the initial gyro-center ion density, *T _{e}* is the initial flux-function equilibrium electron temperature, $K=\mu B+mev\u22a52/2$ is the electron kinetic energy, and $\delta \Phi $ is the potential deviation from the flux-surface averaged part, $\delta \Phi \u2261\Phi \u2212\u27e8\Phi \u27e9$. At the initial time, both $\delta \Phi $ and $\u27e8\Phi \u27e9$ have not developed yet and are zero, and the initial gyrocenter ion density is equal to the real ion density. For $fai$, the Maxwellian distribution without the Boltzmann factor is used.

Since the electron density from $fae$ depends on $\delta \Phi $, the gyrokinetic Poisson equation, Eq. (2) becomes

where the initial flux-function gyrocenter ion density $n\xafi0=\u222bf\xafai$ is equal to *n*_{0} with $f\xaf$ representing the gyro-averaged *f*, $\delta neNA=\u222b(fpe+fge)d3v$ is the electron density that is not represented by the adiabatic Boltzmann distribution function $fae$, and $\delta n\xafi=\u222b(f\xafpi+f\xafgi)d3v$ is the gyro-averaged, perturbed gyrocenter ion density.

To increase the numerical efficiency by reducing data-communication among parallel processors (at the cost of a higher localized arithmetic intensity) in modern massively parallel computers, XGC1 uses an electron sub-cycling scheme that is modified from Ref. 50. The electron time step is smaller by the factor $me/mi$ than the ion time step, and the electric field (thus, the electron charge density) is updated at each ion time step. This sub-cycling scheme is justified by the limitation $\omega \u226ak\u2225ve,th$ that is imposed in the present study. In the original sub-cycling scheme,^{50} the electron marker particle number was much less than the ions: In order to compensate for the statistical error, the electron charge density was evaluated more frequently in inverse-proportion to the number of electrons compared to the ions, and summed up at the ion time step. This required a more frequent data-communication among parallel processors. In the present subcycling scheme, the electron marker particle number is the same as that of ions, and thus the electron charge density is evaluated at the ion time-step.

One issue in the electrostatic simulation of tokamak plasmas is the so-called *ω _{H}* mode.

^{38}The

*ω*mode can give a rapid numerical instability without self-consistent magnetic perturbations unless the time step is small enough to resolve the mode. To avoid this issue, we utilize the numerical scheme of the fluid-kinetic hybrid electron model,

_{H}^{51}utilizing the adiabatic distribution function $fae$ as the fluid part. This scheme has prediction and correction phases. In the prediction phase, the lowest order perturbed potential, $\delta \Phi (0)$ is determined by the adiabatic electron response ( $\delta neNA=0$), and $\delta \Phi (0)$ is used to set $fae$ in the weight evolution equation. In the correction phase, the next order potential $\delta \Phi (1)$ is obtained using the Poisson equation (8) with the non-adiabatic density response $\delta neNA$, which is from the weight evolution equation with $fae(\delta \Phi (0))$. Since the electron particle weights are required only when evaluating $\delta neNA$ for the gyrokinetic Poisson equation, using the direct weight evolution method, the particle weight can be determined according to Eq. (6) only at each ion $\Delta t$ time step when the gyrokinetic Poisson equation is solved.

### B. Linearized Poisson solver, mesh, and charge interpolation

A linearized gyrokinetic Poisson equation is used in this work,

Thus, accuracy for some high amplitude perturbation peaks (XGC1 sees $\delta n/n$ as high as ∼30% in the scrape-off layer), is not guaranteed and their relative error bar can be as high as $O(\delta n/n)2$. Detailed error analysis will be a subject for future study using the recently developed nonlinear Poisson equation solver.

Equation (9) is discretized with a linear finite element method. The flux-surface-average operator in $\delta \Phi =\Phi \u2212\u27e8\Phi \u27e9$ poses a difficulty because it appears as a dense matrix. Thus, Eq. (9) is solved iteratively with $\u27e8\Phi \u27e9$ on the right-hand-side to avoid explicit inversion of the dense flux-surface-average matrix.

An approximately field following mesh (unstructured and triangular) is used. The mesh points are arranged to be axisymmetric to enable more accurate handling of the global axisymmetric physics in the multi-scale simulation of the multiphysics in XGC1. The multiphysics include neoclassical dynamics together with the X-point ion orbit-loss, microturbulence, zonal flows, neutral particle transport, and atomic interaction physics. N-number of “poloidal planes” are defined first at each equidistance in the toroidal angle. Each poloidal plane starts out with the first mesh point on the outboard midplane. Following the magnetic field line in the positive toroidal direction starting from the outboard midplane on the first poloidal plane, we find the intersecting point with the neighboring poloidal plane and mark that as the next mesh point. That mesh point is then mapped to all the poloidal planes axisymmetrically. We continue this “mesh point marking” until the outboard midplane is reached. Unless we are on a mode-rational surface, the last mesh point does not match the first mesh point exactly. We solve this issue by making small adjustments to each mesh point to make the last mesh point to come back to the first mesh point at the outboard midplane. One drawback in this method is that the mesh point distance at the inboard side is denser than that at the outboard side. Near the X-point, the field lines are almost in the toroidal direction ( $B\theta \u21920$). The aforementioned method is not needed since a toroidal-angle following mesh-points is already approximately field-line following.

The toroidal distance ( $2\pi R/N$) between the mesh points is significantly large compared to that on a poloidal plane ( $\u223c\rho i$), where *R* is the major radius and *ρ _{i}* is the ion gyroradius. In order to enhance the accuracy in the charge deposition of marker particles in solving the turbulent fluctuations that are highly stretched in the magnetic field direction ( $k\u2225\u226ak\u22a5$), we execute the charge interpolation in the toroidal direction by following the magnetic field lines on which each marker particle lies.

### C. Verification of the hybrid-Lagrangian scheme with gyro-kinetic ions and drift-kinetic electrons

In Ref. 32, the ion hybrid-Lagrangian scheme was verified using pure ITG turbulence with the adiabatic electron model in a core plasma. The main new feature in the present report is the addition of the kinetic electron hybrid-Lagrangian scheme that is adequate for the simulation of boundary plasma. Before we proceed to the descriptions of the boundary condition for the particle time-marching at the material wall and the dissipative term *S*(*f*), it is prudent at this point to verify the new kinetic electron simulation scheme using the most relevant electrostatic kinetic instability feature in a tokamak plasma; i.e., the transition from the ion-temperature gradient mode to the trapped-electron mode (the so-called “ITG-TEM transition” benchmarking). In the present verification study, we benchmark the XGC1 result against a well-known result by Rewoldt *et al.,*^{52} using the global gyrokinetic codes GTC and GT3D, and a local eigenfrequency code FULL. The gyrokinetic particle code GEM^{53} has participated in this new kinetic-electron benchmarking with XGC1. Figure 1 shows the comparison in the linear growth rate *γ* (left figure) and the real frequency *ω* (right figure) between all five codes. The vertical axis is normalized to $Cs/Lne$, where *C _{s}* is the sound speed and

*L*is the scale-length of the electron density gradient. A good agreement between the XGC1 results and the other code results can be seen. Some deviation of FULL's real frequency in the ITG branch is from the local nature of the eigenfrequency calculation

_{ne}^{52}while other codes are global.

### D. Wall boundary condition for the marker particle time-advance

In XGC1, the wall (divertor and limiter) effects constitute the plasma and solver boundary conditions. There is usually no boundary condition at an inner radius. The plasma particles are absorbed at the wall boundary unless reflected by the sheath potential. All the ions are absorbed, but the electrons whose parallel kinetic energy is lower than the sheath potential are reflected back.

The lost marker particles should still remain in the simulation region (*g* is time invariant) in order to preserve the phase-space volume, but still satisfy *f* = 0 at the wall in the lost particle phase space. In order to achieve these conditions, absorbed marker particles are reflected from the wall with the new weight $w1=\u2212(fa+fg)/(w0g)$ that makes f = 0 at the wall when added to the wall-absorbed particle weight. The full-f weight *w*_{0} is invariant during the time-advancing of the marker particles and is kept invariant here since the reflection process is still part of the *D*/*Dt* time-advancing. The weight change is

We note here that this weight change can create a large magnitude *w*_{1} weight when *f _{g}* = 0 in the beginning of simulation, since

*f*needs to cancel the entire

_{p}*f*upon absorption. However, as the $\Delta w1$ information is gradually transferred to

_{a}*f*in the present hybrid Lagrangian scheme, the effect of the

_{g}*w*

_{1}weight change can be reduced to the fluctuation level.

### E. Logical sheath

In XGC1, the configuration-space grid distance is about half of the ion-gyroradius (∼1 mm). Thus, the sum of Debye sheath, quasi-neutral pre-sheath, and the parallel potential variation that are not resolved by the gyroradius grid needs to be handled as a sub-grid phenomenon: XGC1 uses a modified logical sheath boundary condition.^{55} The total logical sub-grid sheath potential is then determined by counting the number of ions and electrons crossing the wall boundary and by making the lost number of electrons, which has larger parallel energy than the sheath potential, the same as the number of lost ions on average.

Using a thin layer “pill-box” approximation, Gausses law yields the electric field *E _{s}* directed away from the charge density on the wall

*σ*,

where *ϵ*_{0} is the vacuum permittivity. Over a grid spacing $\Delta xg$, this equation becomes

Taking the time derivative and using $d\sigma /dt=e(\Gamma i\u2212\Gamma e)$, where Γ_{i} and Γ_{e} are the absorbed ion and electron fluxes on the wall, we have before further modification $\u2202\Phi s/\u2202t=\u2212(e\Delta xg/\u03f50)(\Gamma i\u2212\Gamma e),$ where $\Phi s$ the logical sheath potential. It is found from numerical trial that this equation brings out a fast oscillation in the value of $\Gamma i\u2212\Gamma e$, including the sign, around the scale of simulation timestep and easily yields a numerical instability.^{56} To avoid this numerical instability, a damping must be imposed on the fast oscillation. We find, empirically, a modified logical sheath algorithm as

where *γ* is the damping factor and *λ _{D}* is the Debye length. We find that $\gamma =10\u22126$ is needed for a stable solution of $\Phi s$. For the initial condition, $e\Phi s=2.5kTe$ is used, which is roughly the Debye sheath potential given by simple analytic theories. Note here that the magnitude of the logical-sheath potential $\Phi s$ is typically greater than the magnitude of the Debye sheath potential (about 1.5 times greater, depending upon the grid size).

### F. Demonstration of edge turbulence simulation across the magnetic separatrix

Before we present the dissipation source terms (nonlinear collisions, heating/cooling, and neutral particles), we demonstrate the capability of the extended hybrid-Lagrangian scheme in simulating electrostatic plasma turbulence across the magnetic separatrix surface. For this demonstration, we turn off the dissipation sources and use only the capabilities and schemes described up to this point.

For this purpose, a DIII-D like L-mode plasma is used with the initial plasma density and temperature profiles modeled as shown in Fig. 2. At the magnetic separatrix, the plasma density is $7\xd71017\u2009m\u22123$ and the electron temperature is 45 eV. This gives 0.003% of electron beta that is smaller than the electron/deuteron mass ratio by about a factor of 10, partially justifying the present electrostatic simulation. A possible electromagnetic effect on edge turbulence will be studied in a future study.

The simulation is performed on the Titan Cray-XK7 supercomputer at the Oak Ridge National Laboratory using 16 384 and some additional reserve compute nodes (90% of Titan capacity) for 72 h. 13 billion marker particles are used for each plasma species. The number of configuration-space grid points of the unstructured triangular mesh is about 55 000 per poloidal plane, with a total of 32 poloidal planes. Because of the high safety factor ( $q95\u22435$) in the edge plasma, 32 poloidal planes is equivalent to $\u2243150$ poloidal planes in the core plasma with q = 1. The velocity space mesh is taken to be a $30\xd731$ rectangular grid with the maximum $v||$ and $v\u22a5$ being three times the thermal velocity. Hence, the total phase space grid has about 1.5 × 10^{9} grid elements. The simulation time step is $1.6\xd710\u22127$ s, and the simulation ended at 1.7 ms in physical time. This is much longer than the saturation time for the collisionless nonlinear turbulence around the magnetic separatrix surface. The simulation time step $1.6\xd710\u22127$ s needs to be small enough to resolve the ion neoclassical orbits, the relevant turbulence modes, and the strongest linear growth time. $1.6\xd710\u22127$ s is about 1/100 of the ion toroidal transit time and less than 1/20 of the maximal relevant turbulence oscillation time scale (300 kHz, see Fig. 8) and the maximal E × B-shearing and linear growth time scale (300 kHz, see Fig. 6). From the convergence study, this time step has been found to be sufficiently small to study the edge turbulence bifurcation physics in the present study.

Figure 3 is the poloidal cross-sectional view of the edge-saturated nonlinear electrostatic turbulence. Core turbulence has a slower growth rate and the growing front can be seen. The inset box shows an enlarged spatial structure of the turbulence. Streamer type structures can be noticed inside the magnetic separatrix (ion temperature gradient and trapped electron turbulence). Around the magnetic separatrix and in the scrape-off layer, coherent blobby type structures can be seen that is intermittent. Detailed statistical properties of this boundary turbulence have been reported elsewhere.^{57} In this subsection, we only present essential figures that have not yet been published, including Figs. 3 and 4.

Figure 4 depicts the square root of the space-time averaged normalized perturbed density $\delta n/n0$ from turbulence, measured at the outboard midplane. Here, *δn* is the perturbation from the toroidally averaged electron density and *n*_{0} is the equilibrium electron density. $(\delta n/n0)2$ is averaged over –30 cm to 30 cm in the poloidal direction, 0 to $2\pi $ in the toroidal direction, and 0.6–1.7 ms in time, when turbulence intensity is nonlinearly saturated. Spatial variation of the flux-surface volume Jacobian is taken into consideration in the averaging operation. The normalized turbulence intensity increases slowly toward the edge, followed by a rapid growth from $\psi n=0.96$, where a large density gradient starts, into the scrape-off layer. The peak of the average turbulence intensity reaches about 9% in the near-scrape-off layer. The individual $\delta n/n0$ fluctuation could reach over 30%, as can be seen from the color bar in Fig. 3. Such a high level turbulence intensity in the boundary plasma, around the separatrix surface into the scrape-off layer, has been observed in multiple tokamaks.^{35}

With the capability demonstration of the hybrid-Lagrangian scheme in the gyrokinetic simulation of the edge turbulence across the magnetic separatrix in this subsection, we will proceed with the description of the additional necessary scheme in Sec. III G: the dissipation terms.

### G. Coulomb collision, heat/cooling source, and neutral atomic physics

The operation *S*(*f*) on the right-hand-side of Eq. (1) includes Coulomb collision, plasma heating, radiative cooling, momentum injection, and neutral atomic physics. The result of the *S*(*f*) operation is reflected on the particle weight *w*_{1} and the grid distribution function *f _{g}* according to the relation $(\Delta w1)w0g+\Delta fg=S(f)\Delta t$ with the time step size $\Delta t$, while the phase variables are unchanged. Note that we hold

*g*unchanged since the particle phase variables are unchanged.

*S*(

*f*) should be operating on the total $f=fa+fg+fp$ and executed on the phase-space grid. The total

*f*for this operation is gathered on the phase space grid by mapping

*f*and

_{p}*f*to the local velocity-space grid. In the simulations performed so far, we choose $(\Delta w1)w0g=S(f)$ when there exist particles on the phase space grid. The modified particle weight is mapped back to the marker particles in the configuration space using the inverse mapping. When there is no particle on the phase-space grid, the change in

_{a}*f*is left on the phase space grid with the modification of

*f*by the amount $\Delta fg=S(f)$.

_{g}The Coulomb collision operator part of *S*(*f*) is solved using a fully nonlinear Fokker-Plank-Landau collision operator on the phase-space grid, i.e., on the 2D velocity space grid at each configuration space grid-node. We refer the readers to Refs. 36 and 48. A fully nonlinear Coulomb collision operation is necessary due to the non-Maxwellian distribution in the edge plasma, especially in the steep pedestal and the scrape-off layer regions, and for the accuracy of the X-point orbit-loss hole physics.

Heating and radiative cooling in *S*(*f*) are applied in an isotropic way in the 2D velocity space. The kinetic energy on each phase-space grid vertex is enhanced or reduced by a small fraction at each simulation time-step, while preserving the particle density and the net parallel momentum on the configuration grid so that the heating and cooling themselves do not inject a net mass or momentum into the plasma. This method preserves the equilibrium (Maxwellian) distribution function. Momentum is injected into the plasma in a similar way, shifting the parallel velocity by a small fraction at each time-step, while preserving the particle density and kinetic energy. The radial profiles of heating, cooling, and momentum sources are either prescribed as input using neutral and impurity profiles. For the cooling profile, an option for a simple radiative power loss model from impurity particles also exists in XGC.

In the XGC version used for the present study, the neutral ionization and charge exchange part of *S*(*f*) is from the built-in Monte Carlo neutral transport routine that uses the rate coefficients of atomic processes.^{58,59} An earlier version XGC0 contains DEGAS2 (Ref. 60) as a subroutine, as described in Ref. 59. First of all, for ions absorbed by the wall, neutral atoms are born with a wall recycling coefficient at the Frank-Condon energy (set to be 3 eV in the present simulation, with random angular distribution) on the material wall or on a specified birth surface in the scrape-off layer (SOL). A gas puff or vacuum pump can also be modeled by generating or annihilating the neutral particles at a prescribed location on the wall.

Again, for the *S*(*f*) operation, all the plasma particles reside on the 5D phase-space grid, while the neutral marker particles reside in the configuration space. As the neutral particles move into the plasma, the charge exchange and ionization probability are calculated using the background plasma information according the atomic rate data, and the new ion marker particles are generated on the ion phase-space grid nodes using the energy and momentum of the neutral particle. The plasma particles are treated as the background fluid (plasma particles are gathered and approximated as shifted Maxwellian on the configuration space grid). Upon the ionization interaction, the new electrons are generated on the electron phase-space grid at the average kinetic energy given in Ref. 61, with a Gaussian velocity distribution. A similar procedure is taken for the plasma marker particle interaction with the background neutral particles. On each phase-space grid point, the charge exchange probability is calculated according the atomic rate data, the new neutral particles are generated using the kinetic energy and momentum on the phase-space grids, and the new ion particle is generated using the neutral particle's energy and momentum. For this plasma atomic process, the neutral particles are treated as background shifted-Maxwellian. The ionization cross-section $\sigma ion$ and the charge-exchange cross-section $\sigma cx$ are as described in^{59}

and the average kinetic energy $Kenew$ of the new electron generated from the electron-impact ionization is as described in^{61}

where $Ke$ is the kinetic energy of the impact electron, $Eionization$ is the ionization energy of the ion species, and the minimum $Kenew$ is set to 4.5 eV. Since the particle-particle interactions are not considered, it is inevitable to suffer from some energy and momentum conservation error in the plasma-neutral system. This topic has been studied in Ref. 59, with the conclusion that the error is not significant enough to replace the built-in neutral Monte-Carlo routine by DEGAS2. The volume recombination has not been considered, assuming that the plasma kinetic energy in the divertor chamber is above the recombination energy.

## IV. GYROKINETIC SIMULATION OF A FAST LOW-TO-HIGH CONFINEMENT MODE BIFURCATION DYNAMICS

### A. Simulation setup and edge plasma profile evolution

A global transport time-scale gyrokinetic investigation of the L-H transition (starting from a global L-mode transport equilibrium, gradually increasing the heating power to get the transition, and observing a pedestal build-up) has been prohibitively expensive on the present-day supercomputers. In the present study, we make the simulation possible by reducing the computational resource requirement as much as possible via a model simplification; i.e., by choosing a fast electrostatic bifurcation case under strong forcing by a high rate of edge heat deposition without prolonging it to the slow, follow-on pedestal build up process. XGC1 simulations and analytic study show that edge turbulence saturation is usually established within 0.1 ms,^{33,62} while in the core plasma, the nonlinear turbulence saturation is established in 1 ms or longer.

By definition, a turbulence-bifurcating plasma is not in a global transport steady state. This implies that the establishment of a global transport steady-state may not be a necessary condition for an edge transport barrier formation study. If the turbulence suppression in the edge layer can occur in much less than 1 ms of plasma time by strong forcing, the transition dynamics can be studied on the 27 peta-flop-peak computer Titan at ORNL.^{63}

For the present study, we use the magnetic field geometry and the plasma profile from the Alcator C-Mod^{34} L-mode plasma discharge #1140613017 as simulation inputs, but taking the toroidal magnetic field (*B _{T}*) to yield the magnetic drift

*V*toward the magnetic X-point (i.e., the favorable direction for an H-mode transition); the actual discharge had

_{B}*V*away from the X-point. The plasma current is parallel to

_{B}*B*. In these plasmas,

_{T}*β*( $\u2261$ the electron kinetic energy density/magnetic energy density) is only 0.01% just inside the separatrix and thus magnetic fluctuation effects are neglected. How the magnetic fluctuations affect the present L-H bifurcation study is a subject of future study, as stated earlier.

_{e}To minimize the computational cost, an exaggerated amount of net heat $\Delta Wlayer\u2248$ 0.8 MW (significantly exceeding the experimentally observed net heat accumulation rate in the edge layer from the 1.6 MW heat flux from core) is accumulated in the $0.947<\psi N<0.989$ edge region so that the edge temperature is forced to increase at an exaggerated rate (Fig. 5), thereby quickly inducing edge transport bifurcation. Here, *ψ _{N}* is a normalized minor radius in terms of poloidal magnetic flux that is zero on the magnetic axis and unity on the separatrix surface. For this plasma, the L-H transition power threshold is known to be about 1–1.5 MW if the magnetic drift is toward the X-point. Thus, the amount of power accumulation in the edge layer is as high as the L-H transition threshold power in this model.

The heat source is designed in the way not to generate an artificial flow in the plasma: After each heating time step in which a small fraction of the particle kinetic energy is increased, any momentum generation is removed by shifting back the particle distribution function in the parallel direction by a proper amount. Moreover, we applied the heat source only at $\psi N<0.76$ so that the heat accumulation in the edge region is from divergence of the radial heat flux.

The edge particle density profile evolves differently from that of temperatures, as can be seen from the bottom figure in Fig. 5. Around the magnetic separatrix surface ( $\psi N>0.98$), a strong neutral ionization raises the particle density, making most of the density slope to exist outside of the separatrix surface. This is caused by the low initial plasma density assumed in the far-scrape-off layer in the simulation, which allows stronger penetration of the neutral particles to the separatrix area. A possible effect of this density buildup around the separatrix on the L-H bifurcation dynamics is unknown in the present study. The rise of particle density deeper in the core at $\psi N<0.95$, while there is no particle source in the core plasma, is due to divergence of turbulent particle flux.

### B. Observation of turbulence bifurcation

A gyrokinetic simulation always experiences transient oscillations in the geodesic acoustic modes (GAMs)^{64–69} as the initial, approximate experimental plasma profile relaxes to a profile that is self-consistent with the gyrokinetic equilibrium and transport. These transient GAMs usually decay away after several oscillations in the near-equilibrium core plasma.^{64–69} However, the GAM oscillation may persist longer or be easily excited in a time-transitional edge plasma due to a weak poloidal winding of the edge magnetic field and a high free energy.^{33,70,71} A strong GAM activity is indeed observed as the L-H bifurcation is approached in ASDEX-U.^{15}

Figure 6(a) depicts the activities of the E × B flow $VE=\u2212Er/B$ (green dashed line), its radial shearing rate $VE\u2032=dVE/dr$ (red dotted), and the turbulence intensity $(\delta n/n)2$ (blue) in the middle of the edge layer at $\psi N\u22480.975$. We will focus our attention to $VE\u2032$, not to *V _{E}* itself, since the latter is found not directly correlated with the bifurcation event. Oscillations at approximately the theoretical GAM frequency of the modeled C-Mod edge plasma ( $\tau GAM\u2243$ 0.03 ms) can be observed in Fig. 6(a). Analysis of these oscillations show that they have an m/n = 0/0 (velocity) and 1/0 (density) Fourier mode structure, consistent with GAMs. The electron and ion heat fluxes are shown in Fig. 6(b) and also exhibit a similar oscillation. Figure 6(c) shows, together with Fig. 6(a), the $VE\u2032$ activities in radius-time with the initial transient shearing rate of $\u223c\u2212105$ Hz decreasing to a negligible level ( $\u223c\u2212104$ Hz) around 0.12 ms at $\psi N\u22430.975$ as the nonlinear edge turbulence is established with its intensity modulated with the global GAM activity. We observe that the GAMs propagate from inner radii towards the edge, with a gradually decreasing radial propagation speed as they approach the edge and some interference pattern as they are reflected from the edge.

A peculiar feature can be noticed in the $VE\u2032$ oscillation before t = 0.175 ms (denoted by the vertical dashed-dotted line): Fig. 6(c) shows that in the edge layer near the magnetic separatrix, the positive peaks of the sheared E × B flow oscillations do not penetrate very well into the region $\psi N>0.96$ (with only one transient penetration, marked with the dashed box), suggesting that at this time there is some mechanism at play to suppress the positive E × B shearing in this region. In the gyrokinetic Poisson equation^{33} $(\rho i2/\lambda i2)VE\u2032\u2248ne\u2212ni,gc$ in an L-mode edge, a negative E × B-flow shearing rate $VE\u2032<0$ implies that the guiding-center plasma is (slightly) positively charged in the edge layer $0.96<\psi N<1$. This also implies that the electrons lead the particle loss, giving rise to the $VE\u2032<0$ before t = 0.175 ms.

Let us continue to discuss the $VE\u2032$ behavior further as seen in Figs. 6(a) and 6(c). At $t\u22480.15$ ms, a strong positive $VE\u2032$ peak penetrates into the far edge region $\psi N\u22430.99$ [Fig. 6(c)] as seen in the dashed box. However, it is only an isolated event. A more robust penetration of $VE\u2032>0$ is observed at $\psi N\u22430.975$. All of a sudden at $t\u22430.175$ ms, $VE\u2032$ becomes positive and increases further in the positive direction [Fig. 6(a)]. $VE\u2032$ and $(\delta n/n)2$ at this position now show an out-of-phase, nonlinear limit cycle behavior. The peak shearing rate at $\psi N\u22480.975$ exceeds 300 kHz at $t\u22480.205$ ms (solid vertical line), which coincides with the maximum linear growth rate of the most unstable dissipative modes^{72} (i.e., dissipative trapped electron modes in the modeled plasma). Also, the second kick into the positive $VE\u2032$ direction that peaks at $t\u22480.205$ ms [Fig. 6(a)] penetrates deeper toward the separatrix $\psi N>0.97$ [Fig. 6(c)]. Around this time, the GAM oscillations at $\psi N<0.95$ are dying out: Thus for $t\u223c0.205\u2009ms$, the stronger penetration of the positive $VE\u2032$ in the region $\psi N>0.97$ is not driven by the stronger GAM activities from the core region. It can also be seen that the sign of the average $VE\u2032$ inside the edge layer $0.96<\psi N<1$ stays positive from $t\u22480.175$ ms, indicating that an electron-dominated particle loss has changed into an ion-dominated loss. Notice here also that the persistent E × B shearing actions are confined to a thin edge layer around $0.96<\psi N<1$.

After $t\u22480.205$ ms, the $VE\u2032$ oscillations at $\psi N\u223c0.975$ cease but $VE\u2032$ grows continuously in the positive direction, and the turbulence is continuously decaying after $\u223c0.22$ ms. The radial thermal fluxes [see Fig. 6(b)] then decay in the same fashion. The ion heat flux does not decay to a small electron heat-flux level because of the significant neoclassical heat flux. The electron particle fluxes also decay similarly as shown in Fig. 7. At this stage, the $VE\u2032>0$ becomes part of the background mean E × B flow shear with a net negative charge (ion-dominated loss). We identify this event as the final stage of turbulence and transport bifurcation after which the pedestal will grow in this fast bifurcation event.

### C. Quenching of turbulence in sheared E × B flow

Theories exit saying that a sheared *E *×* B* flow can (1) shear the turbulence eddies to smaller structures and higher frequency, leading to dissipation at high wave numbers^{73} or (2) quench the turbulence via an eddy tilting-stretching-absorption process.^{9,11} The latter process is different from the former in that it happens via Reynolds work through a conservative process between the turbulence energy and the plasma flow energy. In order to examine how the turbulence behaves in the present simulation during the suppression period, we make contour plots of the turbulence amplitude $\delta \varphi 0$ in the frequency-time space in Fig. 8 at two different radial positions at $\psi N=0.975$ [at the same position where Figs. 6(a) and 6(b) are depicted] and away from the bifurcation layer $\psi N=0.94$.

At $\psi N=0.975$ [corresponding to the white horizontal line in Fig. 6(c)] and $t\u22480.155$ ms (vertical dashed-dotted line), it can be seen that the stronger and lower frequency turbulence begins to decay, and the weaker and higher frequency turbulence begins to appear. The start time of this change coincides with the time when a single strong positive $VE\u2032$ peak of GAM first penetrates to $\psi N\u22480.99$ (marked with a dashed rectangle). The nonlinear limit-cycle behavior starts around $t\u22480.175$ ms [vertical dashed-dotted line, see Fig. 6(a)], and the turbulence intensity suffers a partial cut-off until about 0.19 ms. The observation that the strong transient E × B shearing event at $\psi N\u22480.99$ is felt at $\psi N\u22480.975$ is an indication of a somewhat nonlocal turbulence interaction.

Figure 8(b) shows that a similar phenomena is also observed at some distance away from the edge layer. Even though the main *E *×* B* shearing develops at $\psi N>0.96$, as can be seen in Fig. 6(c), there is a clear sign of turbulence bifurcation activities at $\psi N=0.94$. This phenomenon is not surprising since the radial turbulence wavelength is at least as long as the edge layer width, as shown in Ref. 33 for the edge ITG turbulence and also confirmed from the present simulation for general electrostatic edge turbulence.

In both Figs. 8(a) and 8(b), starting around $t\u22480.205$ ms (solid vertical line), the turbulence begins to get suppressed over the whole frequency range. There is no further cascade of turbulence to any higher frequency than is already reached. This is the time when the E × B shearing rate reaches the critical strength [represented by the solid vertical line in Fig. 6(a)]. The lower frequency turbulence gets suppressed earlier, indicated by the slope starting from the zero-frequency in the time-frequency space. This is a sign that the larger eddies are more susceptible to the E × B shearing effect. These different behaviors before and after 0.205 ms in Fig. 8 suggest that the E × B shearing events suppress the turbulence in two stages.

Since what is observed at $\psi N=0.94$ may better reflect the E × B shearing effect and the turbulence intensity observations presented in Fig. 6(a) due to its proximity, with less interference from the isolated E × B-shearing event at $\psi N\u22430.99$ and $t\u22430.155$ (marked with the dashed box), we attempt to make some detailed comparison between Fig. 8(b) and the time behavior of the Reynolds work rate obtained at $\psi N=0.972$ around the bifurcation time as depicted in Fig. 9. Here, the rate of Reynolds work is plotted as the normalized consumption rate^{7,9,23} $P=\u27e8v\u0303rv\u0303\theta \u27e9VE\u2032/(\gamma effv\u0303\u22a52/2)$ of the turbulence kinetic energy per unit mass, $v\u0303\u22a52/2$, by the *V _{E}* shearing action, using the critical shearing rate 300 kHz as the effective source rate

*γ*of the turbulence kinetic energy (dashed-dotted line). As described earlier in this section, the critical shearing rate has been identified as the strongest growth rate of the turbulence modes (dissipative trapped electron modes in this case). Turbulence intensity from Fig. 6(a) is also plotted for reference (solid line). It is indeed found that the rate of Reynolds work becomes momentarily large enough to consume a significant portion of the turbulence kinetic energy, as indicated by

_{eff}*P*> 1 around $t\u22480.18$ ms and the cut-off of the top in the GAM-oscillating turbulence energy at the corresponding time [in Fig. 6(a)]. Moreover, the time integrated Reynolds work per unit mass ( $5.1\u2009m2/s2$) is somewhat greater than the maximal turbulence energy just before the transition ( $4.5\u2009m2/s2$). We emphasize here that the significant part of the Reynolds work exists only momentarily between the time $t\u22430.175$ ms and $t\u22430.19$ ms. There is a possibility that the turbulence recovers from the driving source, but it does not. From these observations, it can be conjectured that the first-period event is the conservative eddy tilting-stretching-absorption process via Reynolds work, and that the second-period event is the dissipative eddy shearing process from the further growth of the E × B shearing by a mechanism other than the Reynolds work. A more careful study is needed before this conjecture can become a more concrete fact.

There are also differences in the eddy shearing pattern between the two radial locations, as observed from Figs. 8(a) and 8(b). (1) The time slope of the shearing action is lower and the turbulence shearing is saturated at a lower frequency at $\psi N=0.975$ than that at 0.94, and thus (2) the shearing pattern is more coarsely-grained at $\psi N=0.975$ than at 0.94. Despite these differences, the turbulence experiences a similar level of suppression. The difference indicates that the turbulence property at the two places could be different. As depicted in Fig. 10, a visual inspection reveals that the turbulence pattern during the shearing action is coarser around the magnetic separatrix $\psi N>0.97$, where the turbulence pattern appears to be blob-type which is not easily sheared to higher *k _{r}* and frequency, than at inner radii $\psi N<0.96$, where the turbulence is mostly drift-wave type (ITGs-TEMs) with radially elongated streamer structures that can be easily sheared to higher

*k*and frequency. Reference 57 indeed shows from other XGC simulations that the coarser structures around the separatrix and in the scrape-off layer are quasi-coherent and intermittent “blobs.” A figure similar to Fig. 10 that shows the shearing of turbulence from an experimental L-H bifurcation has appeared for a DIII-D discharge.

_{r}^{74}The oppositely directed flows in Fig. 10 across the magnetic separatrix surface are due to a mean E × B flow profile shape resulting from interaction between turbulence, neoclassical (including the ion X-point orbit-loss), neutral particles, and SOL physics at the time of bifurcation.

Figure 11 shows the time-averaged wavenumber spectrum of the turbulence at the two different radial positions before the first-phase E × B shearing starts (t = 0.12–0.17 ms) and well into the second-phase shearing activities (t = 0.22–0.26 ms). The white-red dashed line represents the E × B Doppler-shift $\omega D=2\pi \nu D$ and the white-black dashed line represents the maximum amplitude for each wavenumber. Electron modes are characterized by $(\omega \u2212\omega D)k\theta >0$ and the ion modes are by $(\omega \u2212\omega D)k\theta >0$. Thus, most of the electron modes reside in the positive- $k\theta $ space and the ion modes in the negative- $k\theta $ space for the positive-*ν* values presented here. At t = 0.12–0.17 ms, there are both ion and electron modes at both locations, with the electron mode being stronger at $\psi N=0.941$ than at $\psi N=0.975$. It is interesting to find that the electron modes have already disappeared by t = 0.22–0.26 ms at both locations while the ion modes are being actively sheared away to higher frequency. This may not be a surprising observation since the radially elongated streamers from the trapped electron modes can be highly sensitive to the changes in the drift-resonance of trapped electrons caused by radial E × B shearing. This will be a topic for future research. We note here that the dual modes have been experimentally observed prior to the L-H bifurcation in a low density deuterium discharge in DIII-D.^{74}

### D. Reynolds force and X-point orbit-loss force

Questions that arise at this point include: (1) what triggers the sudden penetration of the strong $VE\u2032>0$ part of the GAM oscillations into the edge layer at $t\u22480.15$ and 0.175 and thus triggering the E × B-shearing of the edge turbulence, (2) why does $VE\u2032$ and its oscillations stay positive and keep on increasing after 0.175 ms, and (3) what maintains the positive E × B flow shear as the turbulence is suppressed?

Figure 6(d), which shows the Reynolds force,^{7,9} $FRe=\u2212d\u27e8v\u0303rv\u0303\theta \u27e9/dr$ at *ψ _{N}* = 0.972 and 0.984, offers an answer to the first question (this force was already hinted by the Reynolds work): there are spatially localized oscillations of

*F*into the positive poloidal direction (electron diamagnetic flow direction) in the edge layer at $t\u22480.155$ and 0.175, with a radial gradient that promotes positive sheared flow in the edge layer (since $dVE/dt\u2248dFRe/dr$). Furthermore, there is another strong peak in the Reynolds force at $t\u22480.195$ at the radial position

_{Re}*ψ*= 0.984, that coincides with another push of the E × B shearing into the positive direction above 300 kHz in Fig. 6(a).

_{N}The second and third questions imply that there is a background force, at least at $t\u22730.175$, pushing the edge layer to a negative charge state or $VE\u2032>0$. The third question also suggests that this background force is strong enough to keep the turbulence suppressed in the edge layer.

The X-point ion orbit loss mechanism^{28,29} is a well-known and robust physics mechanism (see Fig. 12) that drives the edge layer to a negative charge state or a $VE\u2032>0$ state. As the edge *T _{i}* increases, the enlarged X-point ion orbit-loss hole provides a stronger background force leading to a negative local charge with $VE\u2032>0$ that keeps the plasma losses ambipolar. The mechanism can also be interpreted as a loss-hole induced $Jr\xd7B$ return-current force on the main ions that drives a poloidal rotation profile until the viscous force balances the driving $Jr\xd7B$ return-current force in the H-mode equilibrium (see Fig. 8 of Ref. 28).

Figure 6(e) shows a simple estimate of the underlying $Jr\xd7B$ return-current force, measured at $\psi N=0.975$, from the collisionless ion-loss hole in the vicinity of the magnetic X-point as function of time while the local *T _{i}* increases from heating [see Fig. 5(a)]. The orbit-loss driven $Jr\xd7B$ return-current force is comparable and adds to the Reynolds force of Fig. 6(d) in the positive direction. The

*V*behavior in Figs. 6(a) and 6(c), and the second and third questions, can thus be understood as arising from the combined effects of the Reynolds and orbit loss effects.

_{E}Right after the fast turbulence bifurcation, the simulation shows that *E _{r}* (or $VE\u2032$) is not yet supported by the pressure gradient. A strong negative

*dE*/

*dr*at $\psi N\u22481$, or a negative

*E*well in the edge transition layer, has not formed yet either. Assuming that the radial fluid acceleration is small (meaning that the edge plasma is MHD stable), the flux-surface averaged radial force balance equation yields

_{r}Figure 13 plots the poloidal plasma flow speed *V _{p}*, the E × B flow speed $VE\xd7B=\u2212Er/B$, and the ion diamagnetic flow speed $\u2207rP/en$. The figure shows that the

*E*is not supported by the pressure gradient at all in the edge bifurcation layer (

_{r}*E*×

*B*flow is even in the same direction as the ion diamagnetic flow in the bifurcation layer), and that the poloidal flow speed

*V*from the above equation is almost the same as the addition between the

_{p}*E*×

*B*and the ion diamagnetic speed. The toroidal flow speed contribution to

*V*during this fast bifurcation is negligible due to the small $Bp/Bt$. The figure also implies that the negative potential structure and the

_{p}*E*-well have not formed yet, only the curvature has formed (see dotted line in Fig. 14). As the edge pressure profile gradually steepens, it will force

_{r}*E*to be negative in the steep pedestal region to form a negative well, and the usual H-mode pedestal structure is expected to emerge. In other words, across the fast bifurcation time, the most noticeable change in the plasma behavior is in the

_{r}*E*×

*B*shearing rate and the poloidal flow velocity. Change in the toroidal flow speed exists, but not at a level to influence the radial force balance relation yet.

The radial profile of the edge electrostatic potential is depicted in Fig. 14 at *t* = 0.130 ms (solid line), which is before the bifurcation activities start, and at 0.260 ms (dotted line), which is at the end of the bifurcation activities. It can be seen that the curvature flipped sign in the bifurcation layer in the way consistent with the earlier discussions. However, the edge potential well has not formed yet, and the potential value is slightly positive (relative to the grounded wall).

## V. CONCLUSION AND DISCUSSION

A special hybrid-Lagrangian scheme has been developed and extended to enable a total-*δf* gyrokinetic simulation of the boundary-plasma across the closed and the open magnetic field-line topology in the electrostatic limit, for gyrokinetic ions and drift kinetic electrons with a fully nonlinear Fokker-Planck collision operation, and Monte-Carlo neutral atomic particles. The scheme allows ambipolar electron and ion loss to the material wall, which results in neutral atomic particle recycling back into the plasma, but still preserving the phase-space volume of the charged particles. The scheme has been parallelized and optimized to reduce data communications in extreme-scale parallel computers, and applied to study the present L-H bifurcation study,^{31} the divertor heat-flux footprint study,^{62} and other boundary physics.

A fast, forced bifurcation of turbulence and transport has been observed for the first time in an electrostatic nonlinear gyrokinetic simulation. The simulation shows the validity of most of the underlying assumptions used by the popular predator-prey model in this fast transition event, with one important addition that the neoclassical orbit loss physics also plays a critical role in the bifurcation process. We observe that an edge turbulence and transport bifurcation event occurs when the microscale turbulence-driven Reynolds force and the macroscale neoclassical orbit-loss force reinforce each other, and when the combined E × B shearing rate in the edge layer reaches a critical level. Thus, the experimental argument based upon the orbit loss mechanism in Ref. 27 and the conventional Reynolds stress argument work together. In the present simulation of a fast bifurcation event, the Reynolds work causes a conservative eddy tilting-stretching-absorption process, with the orbit-loss then taking over with the shearing of eddies to smaller-scale dissipative eddies and finishing up the turbulence quenching process and keeping the turbulence suppressed.

The present study indicates that an intrinsic limitation of the notion of Reynolds stress in the L-H bifurcation dynamics is its disappearance during the period of turbulence suppression, implying the necessity of some other mechanism for the generation of the sheared E × B flow to keep the turbulence suppressed while a high enough pressure pedestal is formed to provide the needed steady sheared E × B flow. Another limitation is in the lack of a preferred direction in the Reynolds force [it fluctuates in both directions, see Fig. 6(d)]. The synergistic orbit-loss driven E × B-shearing, caused by the rising edge *T _{i}*, that acts in the same direction as the steep $\u2207pi$ driven E × B shearing that develops at a later time, provides such a mechanism. The synergism may help reconcile some experimental observations that ascribe the transition to the orbit loss effect

^{27}or the neoclassical effect

^{30}with reports of the key role of turbulent stress.

^{12–22}There exist other experimental observations that identified a strong correlation between the L-H transition and the orbit loss driven E × B shearing rate.

^{75,76}

The spatial scale of orbit-loss physics is about the ion poloidal gyroradius ( $\Delta \psi N\u22480.05$), while that of the Reynolds stress variation is about $\Delta \psi N\u22480.01$ [see Fig. 6(d)]. The temporal scale of the orbit-loss force development is $\u22480.05$ ms and increasing [Fig. 6(e)], while that of the Reynolds stress is $\u22480.01$ ms [Fig. 6(d)] and fluctuating. Thus, the ion orbit-loss provides a background force, interacting with the space-time dynamical Reynolds force. The ion $90o$ collision time *ν _{ic}* in the transition layer is $\u2248$ 0.05 ms and similar to the ion orbit loss force time scale, and longer than the Reynolds stress time scale. The

*ν*time appears to be related to the limit-cycle time scale [see Fig. 6(a) during $t\u22480.17\u20130.21$ ms].

_{ic}The synergism between the Reynolds and orbit loss forces is also consistent with the general experimental observation that the L-H bifurcation is more difficult in the case when the $\u2207$ B-drift is away from the single-null magnetic X-point, in which the orbit-loss effect is weaker.^{28,29}

In the present simulation, a significant amount of heat ( $\u22480.8$ MW) accumulates in the edge layer in a plasma whose L-H power threshold from the core to the edge is 1–1.5 MW. Thus, the rise of the ion temperature and the orbit-loss force is quick, easily forcing the neoclassical E × B shearing to be effective. In a plasma where the heat deposition to the edge layer is not as strong as what is used in the present simulation, the limit cycle oscillation period may become longer to give a chance for the pedestal pressure gradient to build up, as experimentally observed in Ref. 30. In this case, the E × B shearing from the pressure gradient may add to that from the orbit-loss effect. This is a subject of future study.

## ACKNOWLEDGMENTS

The authors thank L. Delgado-Aparicio for estimating the radiative power loss profile. The authors also thank A. Hubbard for the helpful discussions, especially on the L-I transition. This work was supported mostly by DOE FES and ASCR through the SciDAC project “Partnership Center for Edge Plasma Simulation.” This work was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466, and used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which was supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725 through the ALCC (ASCR Leadership Computing Challenge) program administered by the U.S. DOE Office of Science.