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 ψN0.960.98, where ψN is the normalized poloidal flux, with the time scale 0.1 ms.

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×B-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” model7 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 transition8 and a sharp transition9–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 transition15–20 when operating close to the H-mode power threshold, and a sharp bifurcation21–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 physics28,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 p-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/mi0.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-Mod34 L-mode plasma discharge #1140613017, the plasma βe in the relevant edge region is only 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.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 (60 for deuterons), with ma 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” instability38 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.

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, EFIT40 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.

The electrostatic version XGC1 solves the 5D gyrokinetic Boltzmann equation41,42 in the following form:

ft+Ẋ·fX+v||̇·fv||=S(f)Ẋ=1G[v||b̂+mv||2qB2×b̂+1qB2B×(μBqE¯)]v||̇=1mG(b̂+mv||qB×b̂)·(μBqE¯)G=1+mv||qBb̂·(×b̂).
(1)

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,b̂=B/B, μ=mv2/2B is the magnetic moment, v is the perpendicular velocity to the local magnetic field vector, E¯ 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 

·nemqB2Φ=(1+ρi22)(n¯ine),
(2)

where n¯i is the ion gyro-center density (not the real ion density) averaged over the gyro-orbit, ne is the electron density, is the gradient operator perpendicular to the magnetic field vector B, and ρi is thermal ion gyroradius. For the sake of the discussion in the present paper, we use hydrogenic atoms and replace 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¯ 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

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

DfDt=S(f).
(3)

Note that D/Dt operation itself must conserve the phase space volume42 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+δf, where δf=fp is described by marker particles and f0=fa+fg, with fa being an analytic function and fg residing in the 5D grid-space. Equation (3) can then be equivalently written in the form

DfpDt=Df0Dt+S(f).
(4)

If f0 is set to zero, XGC1 becomes a pure “full-f” particle-in-cell code. If fa is set to be a time-invariant, flux-function Maxwellian fM with fg = 0 and the operator Df0/Dt on the right hand side is simplified to include only the fluctuation-driven term (by neglecting the magnetic B and curvature-drift driven terms), XGC1 becomes a so-called “reduced δf” code with δf=fp. In the reduced δf formalism, a steep radial gradient in fM 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 δffM would be required. These simplifications cannot be enforced in the edge plasma where the magnitude of δ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 δf=fp grows, a small fraction of it can be transferred to fg at every time step time to mitigate the usual “growing-weight problem,”45 as demonstrated in Ref. 32. The Maxwellian part of fg can also be transferred to fM while preserving f0. 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 fp, we use the two-weight scheme.32,46,47 For simplicity, we choose the initial f to be fa, i.e., fg(t=0) =  fp(t=0)=0. The full-f weight w0 is defined as w0gf0(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 w0 defines the density gradient. w0 is invariant in the dissipationless time-advance of the marker particles, i.e., the operation on the left-hand-side of Eq. (4). fp is then the time-advancing marker particle distribution function with the time-advancing δf weight w1, 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

dw1dt=1w0g(Df0Dt+S(f)).
(5)

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

Δw1=Δf0w0g+1wogS(f)dtΔf0w0g+Δtw0gS(f),
(6)

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 f0. This is important especially when we use a noisy interpolated quantity fg 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 fa and fp to the velocity grid and adding to fg. S(f) of Eq. (4) is then evaluated on the 5D continuum grid and its effect on fp 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 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.

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.

The wave frequency ω of interest in XGC1 is much lower than kve,th, where k 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 

fae=n0Te3/2exp(KTe+eδΦTe),
(7)

where n0 is the initial flux-function equilibrium electron density that is equal to the initial gyro-center ion density, Te is the initial flux-function equilibrium electron temperature, K=μB+mev2/2 is the electron kinetic energy, and δΦ is the potential deviation from the flux-surface averaged part, δΦΦΦ. At the initial time, both δΦ and Φ 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 δΦ, the gyrokinetic Poisson equation, Eq. (2) becomes

·nimeB2Φ+(1+ρi22)n0[1exp(δΦTe)]=(1+ρi22)(δn¯iδneNA),
(8)

where the initial flux-function gyrocenter ion density n¯i0=f¯ai is equal to n0 with f¯ representing the gyro-averaged f, δneNA=(fpe+fge)d3v is the electron density that is not represented by the adiabatic Boltzmann distribution function fae, and δn¯i=(f¯pi+f¯gi)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 ωkve,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 ωH 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,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, δΦ(0) is determined by the adiabatic electron response (δneNA=0), and δΦ(0) is used to set fae in the weight evolution equation. In the correction phase, the next order potential δΦ(1) is obtained using the Poisson equation (8) with the non-adiabatic density response δneNA, which is from the weight evolution equation with fae(δΦ(0)). Since the electron particle weights are required only when evaluating δ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 Δt time step when the gyrokinetic Poisson equation is solved.

A linearized gyrokinetic Poisson equation is used in this work,

·n0meB2Φ+(1+ρi22)n0(δΦTe)=(1+ρi22)(δn¯iδneNA).
(9)

Thus, accuracy for some high amplitude perturbation peaks (XGC1 sees δ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(δ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 δΦ=ΦΦ poses a difficulty because it appears as a dense matrix. Thus, Eq. (9) is solved iteratively with Φ 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θ0). The aforementioned method is not needed since a toroidal-angle following mesh-points is already approximately field-line following.

The toroidal distance (2πR/N) between the mesh points is significantly large compared to that on a poloidal plane (ρ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 (kk), we execute the charge interpolation in the toroidal direction by following the magnetic field lines on which each marker particle lies.

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 GEM53 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 Cs is the sound speed and Lne 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 calculation52 while other codes are global.

FIG. 1.

Cross-verification of XGC's new simulation scheme against the known code results for the ITG-TEM mode transition: (a) linear growth rate and (b) real frequency. A reasonably good agreement of the XGC1 results (black X mark) can be seen in comparison with GTC (blue line), GT3D (red line), GEM (purple circle), and the eigenvalue solver Full (black line). The XGC1 and GEM results are overlapped on the FULL, GTC, and GT3D results from Ref. 52. The XGC1, GTC, GT3D, and FULL cross-verification results have already appeared in Ref. 54. Part of this figure is reproduced with permission from Phys. Plasmas 20, 032309 (2013). Copyright 2013 AIP Publishing LLC.

FIG. 1.

Cross-verification of XGC's new simulation scheme against the known code results for the ITG-TEM mode transition: (a) linear growth rate and (b) real frequency. A reasonably good agreement of the XGC1 results (black X mark) can be seen in comparison with GTC (blue line), GT3D (red line), GEM (purple circle), and the eigenvalue solver Full (black line). The XGC1 and GEM results are overlapped on the FULL, GTC, and GT3D results from Ref. 52. The XGC1, GTC, GT3D, and FULL cross-verification results have already appeared in Ref. 54. Part of this figure is reproduced with permission from Phys. Plasmas 20, 032309 (2013). Copyright 2013 AIP Publishing LLC.

Close modal

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=(fa+fg)/(w0g) that makes f = 0 at the wall when added to the wall-absorbed particle weight. The full-f weight w0 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

Δw1=(fa+fg)/(w0g)w1.
(10)

We note here that this weight change can create a large magnitude w1 weight when fg = 0 in the beginning of simulation, since fp needs to cancel the entire fa upon absorption. However, as the Δw1 information is gradually transferred to fg in the present hybrid Lagrangian scheme, the effect of the w1 weight change can be reduced to the fluctuation level.

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 Es directed away from the charge density on the wall σ,

Es=σ/ϵ0,

where ϵ0 is the vacuum permittivity. Over a grid spacing Δxg, this equation becomes

Φs=Δxgσ/ϵ0.

Taking the time derivative and using dσ/dt=e(ΓiΓe), where Γi and Γe are the absorbed ion and electron fluxes on the wall, we have before further modification Φs/t=(eΔxg/ϵ0)(ΓiΓe), where Φs the logical sheath potential. It is found from numerical trial that this equation brings out a fast oscillation in the value of ΓiΓ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

Φst=γeλDϵ0(ΓiΓe),

where γ is the damping factor and λD is the Debye length. We find that γ=106 is needed for a stable solution of Φs. For the initial condition, eΦ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 Φs is typically greater than the magnitude of the Debye sheath potential (about 1.5 times greater, depending upon the grid size).

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×1017m3 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.

FIG. 2.

The modeled initial plasma density, electron temperature, and ion temperature profiles for the capability demonstration of the extended hybrid-Lagrangian scheme in the simulation of the plasma turbulence across the magnetic separatrix surface in a DIII-D like tokamak geometry.

FIG. 2.

The modeled initial plasma density, electron temperature, and ion temperature profiles for the capability demonstration of the extended hybrid-Lagrangian scheme in the simulation of the plasma turbulence across the magnetic separatrix surface in a DIII-D like tokamak geometry.

Close modal

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 (q955) in the edge plasma, 32 poloidal planes is equivalent to 150 poloidal planes in the core plasma with q = 1. The velocity space mesh is taken to be a 30×31 rectangular grid with the maximum v|| and v being three times the thermal velocity. Hence, the total phase space grid has about 1.5 × 109 grid elements. The simulation time step is 1.6×107 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×107 s needs to be small enough to resolve the ion neoclassical orbits, the relevant turbulence modes, and the strongest linear growth time. 1.6×107 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.

FIG. 3.

Poloidal cross-sectional view of the perturbed density normalized to the mean density from electrostatic turbulence in the global boundary plasma. The sheared streamer structures inside the magnetic separatrix surface (black line) and the blobby structures around and outside the magnetic separatrix can be observed. A DIII-D like plasma model is used.

FIG. 3.

Poloidal cross-sectional view of the perturbed density normalized to the mean density from electrostatic turbulence in the global boundary plasma. The sheared streamer structures inside the magnetic separatrix surface (black line) and the blobby structures around and outside the magnetic separatrix can be observed. A DIII-D like plasma model is used.

Close modal
FIG. 4.

Square root of the space-time averaged normalized perturbed density δn/n0 from turbulence, measured at the outboard midplane, where n0 is the equilibrium electron density. (δn/n0)2 is averaged over –30 cm to 30 cm in the poloidal direction, 0 to 2π in the toroidal direction, and 0.6–1.7 ms in time, when turbulence intensity is nonlinearly saturated.

FIG. 4.

Square root of the space-time averaged normalized perturbed density δn/n0 from turbulence, measured at the outboard midplane, where n0 is the equilibrium electron density. (δn/n0)2 is averaged over –30 cm to 30 cm in the poloidal direction, 0 to 2π in the toroidal direction, and 0.6–1.7 ms in time, when turbulence intensity is nonlinearly saturated.

Close modal

Figure 4 depicts the square root of the space-time averaged normalized perturbed density δn/n0 from turbulence, measured at the outboard midplane. Here, δn is the perturbation from the toroidally averaged electron density and n0 is the equilibrium electron density. (δn/n0)2 is averaged over –30 cm to 30 cm in the poloidal direction, 0 to 2π 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 ψ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 δ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.

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 w1 and the grid distribution function fg according to the relation (Δw1)w0g+Δfg=S(f)Δt with the time step size Δ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 fp and fa to the local velocity-space grid. In the simulations performed so far, we choose (Δ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 f is left on the phase space grid with the modification of fg by the amount Δfg=S(f).

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 σion and the charge-exchange cross-section σcx are as described in59 

σion=8×1015Teexp(13.56/Te)1+0.01Tem3s1,σcx=1.1×1014Ei0.3/mim3s1,

and the average kinetic energy Kenew of the new electron generated from the electron-impact ionization is as described in61 

Kenew=Eionization/4,ifKe>1.5Eionization,Kenew=0.5(KeEionization),ifKe<1.5Eionization,

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.

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-Mod34 L-mode plasma discharge #1140613017 as simulation inputs, but taking the toroidal magnetic field (BT) to yield the magnetic drift VB toward the magnetic X-point (i.e., the favorable direction for an H-mode transition); the actual discharge had VB away from the X-point. The plasma current is parallel to BT. In these plasmas, βe ( 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.

To minimize the computational cost, an exaggerated amount of net heat ΔWlayer 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<ψ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.

FIG. 5.

Initial radial profiles of Ti, Te, and ne in the edge region (solid line). The profiles around (t = 0.174 ms, dotted line) and after (t = 0.27 ms, dashed line) the bifurcation event have been plotted together.

FIG. 5.

Initial radial profiles of Ti, Te, and ne in the edge region (solid line). The profiles around (t = 0.174 ms, dotted line) and after (t = 0.27 ms, dashed line) the bifurcation event have been plotted together.

Close modal

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 ψ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 (ψ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 ψN<0.95, while there is no particle source in the core plasma, is due to divergence of turbulent particle flux.

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=Er/B (green dashed line), its radial shearing rate VE=dVE/dr (red dotted), and the turbulence intensity (δn/n)2 (blue) in the middle of the edge layer at ψN0.975. We will focus our attention to VE, not to VE 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 (τGAM 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 activities in radius-time with the initial transient shearing rate of 105 Hz decreasing to a negligible level (104 Hz) around 0.12 ms at ψN0.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.

FIG. 6.

Time behavior at ψN = 0.975 of (a) (δn/n)2, E × B shearing rate, and E × B flow; (b) electron heat flux; (c) E × B shearing rate in radius; (d) Reynolds force at ψN = 0.972 and 0.984; and (e) orbit-loss force. All the one-dimensional quantities are flux surface averaged values except (δn/n)2 in (a), which is measured at the outboard midplane. This figure is reproduced with permission from Chang et al., Phys. Rev. Lett. 118, 175001 (2017). Copyright 2017 American Physical Society.

FIG. 6.

Time behavior at ψN = 0.975 of (a) (δn/n)2, E × B shearing rate, and E × B flow; (b) electron heat flux; (c) E × B shearing rate in radius; (d) Reynolds force at ψN = 0.972 and 0.984; and (e) orbit-loss force. All the one-dimensional quantities are flux surface averaged values except (δn/n)2 in (a), which is measured at the outboard midplane. This figure is reproduced with permission from Chang et al., Phys. Rev. Lett. 118, 175001 (2017). Copyright 2017 American Physical Society.

Close modal

A peculiar feature can be noticed in the VE 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 ψ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 equation33 (ρi2/λi2)VEneni,gc in an L-mode edge, a negative E × B-flow shearing rate VE<0 implies that the guiding-center plasma is (slightly) positively charged in the edge layer 0.96<ψN<1. This also implies that the electrons lead the particle loss, giving rise to the VE<0 before t = 0.175 ms.

Let us continue to discuss the VE behavior further as seen in Figs. 6(a) and 6(c). At t0.15 ms, a strong positive VE peak penetrates into the far edge region ψN0.99 [Fig. 6(c)] as seen in the dashed box. However, it is only an isolated event. A more robust penetration of VE>0 is observed at ψN0.975. All of a sudden at t0.175 ms, VE becomes positive and increases further in the positive direction [Fig. 6(a)]. VE and (δn/n)2 at this position now show an out-of-phase, nonlinear limit cycle behavior. The peak shearing rate at ψN0.975 exceeds 300 kHz at t0.205 ms (solid vertical line), which coincides with the maximum linear growth rate of the most unstable dissipative modes72 (i.e., dissipative trapped electron modes in the modeled plasma). Also, the second kick into the positive VE direction that peaks at t0.205 ms [Fig. 6(a)] penetrates deeper toward the separatrix ψN>0.97 [Fig. 6(c)]. Around this time, the GAM oscillations at ψN<0.95 are dying out: Thus for t0.205ms, the stronger penetration of the positive VE in the region ψ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 inside the edge layer 0.96<ψN<1 stays positive from t0.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<ψN<1.

After t0.205 ms, the VE oscillations at ψN0.975 cease but VE grows continuously in the positive direction, and the turbulence is continuously decaying after 0.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>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.

FIG. 7.

Flux-surface-averaged radial particle flux at ψN=0.975 corresponding to Fig. 6(b), showing the suppression around t = 0.175 ms.

FIG. 7.

Flux-surface-averaged radial particle flux at ψN=0.975 corresponding to Fig. 6(b), showing the suppression around t = 0.175 ms.

Close modal

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 numbers73 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 δϕ0 in the frequency-time space in Fig. 8 at two different radial positions at ψN=0.975 [at the same position where Figs. 6(a) and 6(b) are depicted] and away from the bifurcation layer ψN=0.94.

FIG. 8.

Contour plot of the turbulence intensity ϕ at the outboard midplane in the frequency-time space at two different radial positions: (a) at ψN=0.975 [at the same position where Figs. 6(a) and 6(b) are depicted] and (b) away from the bifurcation layer ψN=0.94.

FIG. 8.

Contour plot of the turbulence intensity ϕ at the outboard midplane in the frequency-time space at two different radial positions: (a) at ψN=0.975 [at the same position where Figs. 6(a) and 6(b) are depicted] and (b) away from the bifurcation layer ψN=0.94.

Close modal

At ψN=0.975 [corresponding to the white horizontal line in Fig. 6(c)] and t0.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 peak of GAM first penetrates to ψN0.99 (marked with a dashed rectangle). The nonlinear limit-cycle behavior starts around t0.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 ψN0.99 is felt at ψN0.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 ψN>0.96, as can be seen in Fig. 6(c), there is a clear sign of turbulence bifurcation activities at ψ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 t0.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 ψ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 ψN0.99 and t0.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 ψN=0.972 around the bifurcation time as depicted in Fig. 9. Here, the rate of Reynolds work is plotted as the normalized consumption rate7,9,23 P=ṽrṽθVE/(γeffṽ2/2) of the turbulence kinetic energy per unit mass, ṽ2/2, by the VE shearing action, using the critical shearing rate 300 kHz as the effective source rate γeff 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 P > 1 around t0.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.1m2/s2) is somewhat greater than the maximal turbulence energy just before the transition (4.5m2/s2). We emphasize here that the significant part of the Reynolds work exists only momentarily between the time t0.175 ms and t0.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.

FIG. 9.

Normalized consumption rate ṽrṽθV/(γeffṽ2/2) of the turbulence kinetic energy ṽ2/2 by the shearing action VE at ψN=0.972 around the bifurcation time, using the critical shearing rate 300 kHz as the effective source rate (γeff). Turbulence intensity dynamics is also shown. These quantities are measured at the outboard midplane. This figure is reproduced with permission from Chang et al., Phys. Rev. Lett. 118, 175001 (2017). Copyright 2017 American Physical Society.

FIG. 9.

Normalized consumption rate ṽrṽθV/(γeffṽ2/2) of the turbulence kinetic energy ṽ2/2 by the shearing action VE at ψN=0.972 around the bifurcation time, using the critical shearing rate 300 kHz as the effective source rate (γeff). Turbulence intensity dynamics is also shown. These quantities are measured at the outboard midplane. This figure is reproduced with permission from Chang et al., Phys. Rev. Lett. 118, 175001 (2017). Copyright 2017 American Physical Society.

Close modal

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 ψN=0.975 than that at 0.94, and thus (2) the shearing pattern is more coarsely-grained at ψ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 ψN>0.97, where the turbulence pattern appears to be blob-type which is not easily sheared to higher kr and frequency, than at inner radii ψ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 kr 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.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.

FIG. 10.

E × B shearing action of the turbulence eddies across the magnetic separatrix surface.

FIG. 10.

E × B shearing action of the turbulence eddies across the magnetic separatrix surface.

Close modal

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 ωD=2πνD and the white-black dashed line represents the maximum amplitude for each wavenumber. Electron modes are characterized by (ωωD)kθ>0 and the ion modes are by (ωωD)kθ>0. Thus, most of the electron modes reside in the positive- kθ space and the ion modes in the negative- kθ 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 ψN=0.941 than at ψ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 

FIG. 11.

Time-averaged wavenumber spectrum of the turbulence at the two different radial positions (ψN=0.941 and 0.975) measured at the outboard midplane 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 ωD=2πνD and the white-black dashed line represents the maximum amplitude for each wavenumber. Electron modes are characterized by (ωωD)kθ>0 and the ion modes are by (ωωD)kθ>0. Thus, most of the electron modes reside in the positive- kθ space and the ion modes in the negative- kθ space for positive-ν values are presented here.

FIG. 11.

Time-averaged wavenumber spectrum of the turbulence at the two different radial positions (ψN=0.941 and 0.975) measured at the outboard midplane 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 ωD=2πνD and the white-black dashed line represents the maximum amplitude for each wavenumber. Electron modes are characterized by (ωωD)kθ>0 and the ion modes are by (ωωD)kθ>0. Thus, most of the electron modes reside in the positive- kθ space and the ion modes in the negative- kθ space for positive-ν values are presented here.

Close modal

Questions that arise at this point include: (1) what triggers the sudden penetration of the strong VE>0 part of the GAM oscillations into the edge layer at t0.15 and 0.175 and thus triggering the E × B-shearing of the edge turbulence, (2) why does VE 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=dṽrṽθ/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 FRe into the positive poloidal direction (electron diamagnetic flow direction) in the edge layer at t0.155 and 0.175, with a radial gradient that promotes positive sheared flow in the edge layer (since dVE/dtdFRe/dr). Furthermore, there is another strong peak in the Reynolds force at t0.195 at the radial position ψN = 0.984, that coincides with another push of the E × B shearing into the positive direction above 300 kHz in Fig. 6(a).

The second and third questions imply that there is a background force, at least at t0.175, pushing the edge layer to a negative charge state or VE>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 mechanism28,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>0 state. As the edge Ti increases, the enlarged X-point ion orbit-loss hole provides a stronger background force leading to a negative local charge with VE>0 that keeps the plasma losses ambipolar. The mechanism can also be interpreted as a loss-hole induced Jr×B return-current force on the main ions that drives a poloidal rotation profile until the viscous force balances the driving Jr×B return-current force in the H-mode equilibrium (see Fig. 8 of Ref. 28).

FIG. 12.

Two types of X-loss orbits. This figure is reproduced with permission from Phys. Plasmas 11, 5626 (2004). Copyright 2004 AIP Publishing LLC.77 

FIG. 12.

Two types of X-loss orbits. This figure is reproduced with permission from Phys. Plasmas 11, 5626 (2004). Copyright 2004 AIP Publishing LLC.77 

Close modal

Figure 6(e) shows a simple estimate of the underlying Jr×B return-current force, measured at ψN=0.975, from the collisionless ion-loss hole in the vicinity of the magnetic X-point as function of time while the local Ti increases from heating [see Fig. 5(a)]. The orbit-loss driven Jr×B return-current force is comparable and adds to the Reynolds force of Fig. 6(d) in the positive direction. The VE 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.

Right after the fast turbulence bifurcation, the simulation shows that Er (or VE) is not yet supported by the pressure gradient. A strong negative dE/dr at ψN1, or a negative Er 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

Vp=rP/enEr/B+VtBp/Bt.

Figure 13 plots the poloidal plasma flow speed Vp, the E × B flow speed VE×B=Er/B, and the ion diamagnetic flow speed rP/en. The figure shows that the Er is not supported by the pressure gradient at all in the edge bifurcation layer (E × B flow is even in the same direction as the ion diamagnetic flow in the bifurcation layer), and that the poloidal flow speed Vp from the above equation is almost the same as the addition between the E × B and the ion diamagnetic speed. The toroidal flow speed contribution to Vp during this fast bifurcation is negligible due to the small Bp/Bt. The figure also implies that the negative potential structure and the Er-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 Er 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 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.

FIG. 13.

Flux-surface-averaged poloidal plasma flow, E × B flow, and the ion diamagnetic flow in the edge layer, showing that the Er is not generated by the pressure gradient during this fast bifurcation event.

FIG. 13.

Flux-surface-averaged poloidal plasma flow, E × B flow, and the ion diamagnetic flow in the edge layer, showing that the Er is not generated by the pressure gradient during this fast bifurcation event.

Close modal
FIG. 14.

Electrostatic potential at the outboard midplane as function of normalized radius ψN before the bifurcation activities start (t = 0.130 ms) and near the end of the bifurcation activities (t = 0.260 ms).

FIG. 14.

Electrostatic potential at the outboard midplane as function of normalized radius ψN before the bifurcation activities start (t = 0.130 ms) and near the end of the bifurcation activities (t = 0.260 ms).

Close modal

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).

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 Ti, that acts in the same direction as the steep pi 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 effect27 or the neoclassical effect30 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 (ΔψN0.05), while that of the Reynolds stress variation is about ΔψN0.01 [see Fig. 6(d)]. The temporal scale of the orbit-loss force development is 0.05 ms and increasing [Fig. 6(e)], while that of the Reynolds stress is 0.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 0.05 ms and similar to the ion orbit loss force time scale, and longer than the Reynolds stress time scale. The νic time appears to be related to the limit-cycle time scale [see Fig. 6(a) during t0.170.21 ms].

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 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 (0.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.

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.

1.
A.
Or
and
F.
Busse
,
J. Fluid Mech.
174
,
313
(
1987
).
2.
M.
McIntyre
,
J. Atmos. Terr. Phys.
51
,
29
(
1989
).
3.
F.
Wagner
,
G.
Becker
,
K.
Behringer
,
D.
Campbell
,
A.
Eberhagen
,
W.
Engelhardt
,
G.
Fussmann
,
O.
Gehre
,
J.
Gernhardt
,
G. V.
Gierke
 et al,
Phys. Rev. Lett.
49
,
1408
(
1982
).
4.
C.
Gormezano
,
A.
Sips
,
T.
Luce
,
S.
Ide
,
A.
Becoulet
,
X.
Litaudon
,
A.
Isayama
,
J.
Hobirk
,
M.
Wade
,
T.
Oikawa
 et al,
Nucl. Fusion
47
,
S285
(
2007
).
5.
K.
Ikeda
,
Nucl. Fusion
47
(
2007
).
6.
F.
Wagner
,
Plasma Phys. Controlled Fusion
49
,
B1
(
2007
).
7.
E.
Kim
and
P.
Diamond
,
Phys. Rev. Lett.
90
,
185006
(
2003
).
8.
W.
Weymiens
,
H.
de Blank
,
G.
Hogeweij
, and
J. C.
de Valena
,
Phys. Plasmas
19
,
072309
(
2012
).
9.
P.
Diamond
,
S.-I.
Itoh
,
K.
Itoh
, and
T.
Hahm
,
Plasma Phys. Controlled Fusion
47
,
R35
(
2005
).
10.
K.
Miki
,
P. H.
Diamond
,
D.
Grcan
,
G. R.
Tynan
,
T.
Estrada
,
L.
Schmitz
, and
G. S.
Xu
,
Phys. Plasmas
19
,
092306
(
2012
).
11.
G.
Tynan
,
I.
Cziegler
,
P.
Diamond
,
M.
Malkov
,
A.
Hubbard
,
J.
Hughes
,
J.
Terry
, and
J.
Irby
,
Plasma Phys. Controlled Fusion
58
,
044003
(
2016
).
12.
L.
Chn
,
P.
Beyer
,
Y.
Sarazin
,
G.
Fuhr
,
C.
Bourdelle
, and
S.
Benkadda
,
Phys. Plasmas
21
,
070702
(
2014
).
13.
G. Y.
Park
,
S. S.
Kim
,
H.
Jhang
,
P. H.
Diamond
,
T.
Rhee
, and
X. Q.
Xu
,
Phys. Plasmas
22
,
032505
(
2015
).
14.
J. J.
Rasmussen
,
A.
Nielsen
,
J.
Madsen
,
V.
Naulin
, and
G. S.
Xu
,
Plasma Phys. Controlled Fusion
58
,
014031
(
2016
).
15.
G. D.
Conway
,
C.
Angioni
,
F.
Ryter
,
P.
Sauter
,
J.
Vicente
, and
ASDEX Upgrade Team
,
Phys. Rev. Lett.
106
,
065001
(
2011
).
16.
T.
Estrada
,
T.
Happel
,
C.
Hidalgo
,
E.
Ascasbar
, and
E.
Blanco
,
Europhys. Lett.
92
,
35001
(
2010
).
17.
P.
Manz
,
M.
Ramisch
, and
U.
Stroth
,
Phys. Rev. E
82
,
056403
(
2010
).
18.
G. S.
Xu
,
B. N.
Wan
,
H. Q.
Wang
,
H. Y.
Guo
,
H. L.
Zhao
,
A. D.
Liu
,
V.
Naulin
,
P. H.
Diamond
,
G. R.
Tynan
,
M.
Xu
 et al,
Phys. Rev. Lett.
107
,
125001
(
2011
).
19.
L.
Schmitz
,
L.
Zeng
,
T. L.
Rhodes
,
J. C.
Hillesheim
,
E. J.
Doyle
,
R. J.
Groebner
,
W. A.
Peebles
,
K. H.
Burrell
, and
G.
Wang
,
Phys. Rev. Lett.
108
,
155002
(
2012
).
20.
G.
Tynan
,
M.
Xu
,
P.
Diamond
,
J.
Boedo
,
I.
Cziegler
,
N.
Fedorczak
,
P.
Manz
,
K.
Miki
,
S.
Thakur
,
L.
Schmitz
 et al,
Nucl. Fusion
53
,
073053
(
2013
).
21.
C.
Hidalgo
,
M. A.
Pedrosa
,
E.
Snchez
,
R.
Balbn
,
A.
Lpez-Fraguas
,
B.
van Milligen
,
C.
Silva
,
H.
Fernandes
,
C. A. F.
Varandas
,
C.
Riccardi
 et al,
Plasma Phys. Controlled Fusion
42
,
A153
(
2000
).
22.
G.
McKee
,
P.
Gohil
,
D.
Schlossberg
,
J.
Boedo
,
K.
Burrell
,
J.
deGrassie
,
R.
Groebner
,
R.
Moyer
,
C.
Petty
,
T.
Rhodes
 et al,
Nucl. Fusion
49
,
115016
(
2009
).
23.
P.
Manz
,
G. S.
Xu
,
B. N.
Wan
,
H. Q.
Wang
,
H. Y.
Guo
,
I.
Cziegler
,
N.
Fedorczak
,
C.
Holland
,
S. H.
Mller
,
S. C.
Thakur
 et al,
Phys. Plasmas
19
,
072311
(
2012
).
24.
Z.
Yan
,
G. R.
McKee
,
R.
Fonck
,
P.
Gohil
,
R. J.
Groebner
, and
T. H.
Osborne
,
Phys. Rev. Lett.
112
,
125002
(
2014
).
25.
I.
Cziegler
,
P. H.
Diamond
,
N.
Fedorczak
,
P.
Manz
,
G. R.
Tynan
,
M.
Xu
,
R. M.
Churchill
,
A. E.
Hubbard
,
B.
Lipschultz
,
J. M.
Sierchio
 et al,
Phys. Plasmas
20
,
055904
(
2013
).
26.
R. A.
Moyer
,
K. H.
Burrell
,
T. N.
Carlstrom
,
S.
Coda
,
R. W.
Conn
,
E. J.
Doyle
,
P.
Gohil
,
R. J.
Groebner
,
J.
Kim
,
R.
Lehmer
 et al,
Phys. Plasmas
2
,
2397
(
1995
).
27.
T.
Kobayashi
,
K.
Itoh
,
T.
Ido
,
K.
Kamiya
,
S.-I.
Itoh
,
Y.
Miura
,
Y.
Nagashima
,
A.
Fujisawa
,
S.
Inagaki
,
K.
Ida
 et al,
Phys. Rev. Lett.
111
,
035002
(
2013
).
28.
C. S.
Chang
,
S.
Kue
, and
H.
Weitzner
,
Phys. Plasmas
9
,
3884
(
2002
).
29.
C. S.
Chang
,
S.
Ku
, and
H.
Weitzner
,
Phys. Plasmas
11
,
2649
(
2004
).
30.
M.
Cavedon
,
T.
Ptterich
,
E.
Viezzer
,
G.
Birkenmeier
,
T.
Happel
,
F. M.
Laggner
,
P.
Manz
,
F.
Ryter
,
U.
Stroth
, and
ASDEX Upgrade Team
,
Nucl. Fusion
57
,
014002
(
2017
).
31.
C. S.
Chang
,
S.
Ku
,
G. R.
Tynan
,
R.
Hager
,
R. M.
Churchill
,
I.
Cziegler
,
M.
Greenwald
,
A. E.
Hubbard
, and
J. W.
Hughes
,
Phys. Rev. Lett.
118
,
175001
(
2017
).
32.
S.
Ku
,
R.
Hager
,
C.
Chang
,
J.
Kwon
, and
S.
Parker
,
J. Comput. Phys.
315
,
467
(
2016
).
33.
C. S.
Chang
,
S.
Ku
,
P. H.
Diamond
,
Z.
Lin
,
S.
Parker
,
T. S.
Hahm
, and
N.
Samatova
,
Phys. Plasmas
16
,
056108
(
2009
).
34.
E. S.
Marmar
and
Alcator C-Mod Group
,
Fusion Sci. Technol.
51
,
261
(
2007
).
35.
J. M. D. A.
D'Ippolito
and
S.
Zweben
,
Phys. Plasmas
18
,
060501
(
2011
).
36.
R.
Hager
,
E.
Yoon
,
S.
Ku
,
E.
D'Azevedo
,
P.
Worley
, and
C.
Chang
,
J. Comput. Phys.
315
,
644
(
2016
), ISSN: 0021-9991.
37.
R.
Courant
,
K.
Friedrichs
, and
H.
Lewy
,
Math. Ann.
100
,
32
(
1928
).
38.
W.
Lee
,
J. Comput. Phys.
72
,
243
(
1987
).
39.
S.
Ku
,
C.
Chang
, and
P.
Diamond
,
Nucl. Fusion
49
,
115021
(
2009
).
40.
L.
Lao
,
H. S.
John
,
R.
Stambaugh
,
A.
Kellman
, and
W.
Pfeiffer
,
Nucl. Fusion
25
,
1611
(
1985
).
41.
R. G.
Littlejohn
,
Phys. Fluids
28
,
2015
(
1985
).
42.
T. S.
Hahm
,
Phys. Fluids
31
,
2670
(
1988
).
43.
P. N.
Guzdar
,
C. S.
Liu
,
J. Q.
Dong
, and
Y. C.
Lee
,
Phys. Rev. Lett.
57
,
2818
(
1986
).
44.
Y. C.
Lee
,
J. Q.
Dong
,
P. N.
Guzdar
, and
C. S.
Liu
,
Phys. Fluids
30
,
1331
(
1987
).
45.
J. A.
Krommes
and
G.
Hu
,
Phys. Plasmas
1
,
3211
(
1994
).
46.
Y.
Chen
and
R.
White
,
Phys. Plasmas
4
,
3591
(
1997
).
47.
W.
Wang
,
T.
Hahm
,
W.
Lee
,
G.
Rewoldt
,
J.
Manickam
, and
W.
Tang
,
Phys. Plasmas
14
,
072306
(
2007
).
48.
E.
Yoon
and
C.
Chang
,
Phys. Plasmas
21
,
032503
(
2014
).
49.
I.
Manuilskiy
and
W. W.
Lee
,
Phys. Plasmas
7
,
1381
(
2000
).
50.
J.
Adam
,
A.
Courdin
, and
A.
Langdon
,
J. Comput. Phys.
47
,
229
(
1982
).
51.
Z.
Lin
and
L.
Chen
,
Phys. Plasmas
8
,
1447
(
2001
).
52.
G.
Rewoldt
,
Z.
Lin
, and
Y.
Idomura
,
Comput. Phys. Commun.
177
,
775
(
2007
), ISSN: 0010-4655.
53.
Y.
Chen
and
S.
Parker
,
Phys. Plasmas
8
,
2095
(
2001
).
54.
I.
Holod
and
Z.
Lin
,
Phys. Plasmas
20
,
032309
(
2013
).
55.
S.
Parker
,
R.
Procassini
,
C.
Birdsall
, and
B.
Cohen
,
J. Comput. Phys.
104
,
41
(
1993
).
56.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics via Computer Simulation
(
McGraw-Hill
,
1985
), ISBN: 0-07-005371-5.
57.
R. M.
Churchill
,
C. S.
Chang
,
S.
Ku
, and
J.
Dominski
,
Plasma Phys. Controlled Fusion
59
,
105014
(
2017
).
58.
J.
Huba
, see https://www.nrl.navy.mil/ppd/content/nrl-plasma-formulary for NRL Plasma Formulary (
1994
).
59.
D.
Stotler
,
J.
Lang
,
C.
Chang
,
R.
Churchill
, and
S.
Ku
,
Nucl. Fusion
57
,
086028
(
2017
).
60.
D. P.
Stotler
and
C. F. F.
Karney
,
Contrib. Plasma Phys.
34
,
392
(
1994
).
61.
R. K.
Janev
,
W. D.
Langer
,
J. K.
Evans
, and
J. D. E.
Post
,
Elementary Processes in Hydrogen-Helium Plasmas
, Springer Series on Atoms and Plasmas (
Springer-Verlag
,
New York
,
1987
).
62.
C.
Chang
,
S.
Ku
,
A.
Loarte
,
V.
Parail
,
F.
Kchl
,
M.
Romanelli
,
R.
Maingi
,
J.-W.
Ahn
,
T.
Gray
,
J.
Hughes
 et al,
Nucl. Fusion
57
,
116023
(
2017
).
63.
See https://www.olcf.ornl.gov/computing?resources/titan?cray?xk7/ for Titan supercomputer at Oak Ridge National Laboratory.
64.
M. N.
Rosenbluth
and
F. L.
Hinton
,
Phys. Rev. Lett.
80
,
724
(
1998
).
65.
N.
Winsor
,
J. L.
Johnson
, and
J. M.
Dawson
,
Phys. Fluids
11
,
2448
(
1968
).
66.
S. V.
Novakovskii
,
C. S.
Liu
,
R. Z.
Sagdeev
, and
M. N.
Rosenbluth
,
Phys. Plasmas
4
,
4272
(
1997
).
67.
V. B.
Lebedev
,
P. N.
Yushmanov
,
P. H.
Diamond
,
S. V.
Novakovskii
, and
A. I.
Smolyakov
,
Phys. Plasmas
3
,
3023
(
1996
).
68.
H.
Sugama
and
T.-H.
Watanabe
,
J. Plasma Phys.
72
,
825
(
2006
).
69.
F.
Zonca
and
L.
Chen
,
Europhys. Lett.
83
,
35001
(
2008
).
70.
R.
Hager
and
K.
Hallatschek
,
Phys. Rev. Lett.
108
,
035004
(
2012
).
71.
K.
Hallatschek
and
D.
Biskamp
,
Phys. Rev. Lett.
86
,
1223
(
2001
).
72.
M.
Romanelli
,
G.
Regnoli
, and
C.
Bourdelle
,
Phys. Plasmas
14
,
082305
(
2007
).
73.
H.
Biglari
,
P. H.
Diamond
, and
P.
Terry
,
Phys. Fluids B
2
,
1
(
1990
).
74.
Z.
Yan
,
P.
Gohil
,
G.
McKee
,
D.
Eldon
,
B.
Grierson
,
T.
Rhodes
, and
C.
Petty
,
Nucl. Fusion
57
,
126015
(
2017
).
75.
D.
Battaglia
,
C.
Chang
,
S.
Kaye
,
K.
Kim
,
S.
Ku
,
R.
Maingi
,
R.
Bell
,
A.
Diallo
,
S.
Gerhardt
,
B.
LeBlanc
 et al,
Nucl. Fusion
53
,
113032
(
2013
).
76.
S.
Kaye
,
R.
Maingi
,
D.
Battaglia
,
R.
Bell
,
C.
Chang
,
J.
Hosea
,
H.
Kugel
,
B.
LeBlanc
,
H.
Meyer
,
G.
Park
 et al,
Nucl. Fusion
51
,
113019
(
2011
).
77.
S.
Ku
,
H.
Baek
, and
C. S.
Chang
,
Phys. Plasmas
11
,
5626
(
2004
).