Full-f version of GENE for turbulence in open-field-line systems

Understanding and controlling properties of plasmas in the tokamak edge, including the pedestal and the scrape-off layer (SOL) separated by the last-closed flux surface (LCFS), is crucial for achieving successful magnetic confinement, due to the fact that plasma properties (density and temperature) at the edge significantly influence core confinement, and the high power load in a narrow layer outside of LCFS can jeopardize the plasma-facing-components (PFCs) of the divertor or limiter plates.^1–4 Edge plasmas are complex and present unique features that are typically absent in core plasmas, such as steep temperature and density gradients at the pedestal, large amplitudes of disturbances and intermittent structures (both electrostatic and electromagnetic), plasma–wall interactions, as well as complicated X-point geometry that includes topologically separated open and closed field lines.^5–7 Major efforts are being devoted to developing gyrokinetic particle-in-cell codes with capabilities of simulating the SOL region,^8,9 and new gyrokinetic continuum codes designed for edge plasmas.^10–13 Meanwhile, comprehensive and sophisticated gyrokinetic codes have been developed for simulating core turbulence and explaining relevant experimental data in the past two decades or so.^14–21 To various degrees, the aforementioned characteristics of edge plasmas pose challenges and questions when core turbulence codes are applied to modeling edge plasmas, and there have been few efforts to modify these codes to handle them.

In this work, we extend the global version of the 3D2V core turbulence code GENE,^14,22,23 to model large-amplitude electrostatic disturbances in open-field-line systems. Large-amplitude disturbances and open boundaries are major characteristics of the SOL. This study is meaningful in three aspects. First, the GENE code employs a δf-splitting technique to achieve high numerical accuracy for core turbulence simulations. The parallel nonlinearity is usually neglected due to the classical gyrokinetic ordering $δf/F≪1$ in the tokamak core.²⁴ The potential breakdown of this ordering results in questioning the applicability of δf-splitting in simulating edge plasmas with large-amplitude disturbances.⁶ As shown in our previous study²⁵ based on a 1D1V code extracted from the 3D2V GENE code, by including the parallel nonlinearity and combining it with the corresponding linear term in a finite-volume discretization, the δf-splitting causes no numerical difficulties in handling the large-amplitude disturbances. In this study, the same approach is taken for the 3D2V δf code, which then becomes equivalent to and can be run as a full-f code. Second, an open parallel boundary is implemented to account for plasma sheath effects in gyrokinetic simulations. The open boundary takes a conducting-sheath boundary model, in which the particle loss to and reflection from the wall are determined by the electrostatic potential at the simulation domain edge, and net currents are allowed to flow through the boundary freely. This sheath model was introduced by Shi et al.¹² and implemented in the Gkeyll code, which features a discontinuous Galerkin (finite-element) method. In this study, this model is implemented in the context of a finite-volume method. Since plasma sheath beyond the boundary is not resolved, the distributions and potential are reconstructed by using local polynomials to interpolate the data inside the domain, in order to evaluate the outgoing particle flux and the parallel nonlinearity near the boundary. Third, the Lenard–Bernstein collision operator acting on full-f is implemented. The implementation utilizes a finite-volume discretization structure that has been used for the linearized Fokker-Planck collision operator acting on δf. The number density and energy are conserved numerically. The collision operator contains pitch-angle scattering and velocity-space diffusion that preferentially damps small scale structures. The present extensions provide GENE with essential capabilities to simulate SOL plasmas and constitute an initial effort for future comprehensive modeling of edge plasmas.

To test the extensions, the code is applied to a test case of electrostatic plasma turbulence in the Large Plasma Device (LAPD)²⁶ at UCLA. The LAPD plasma contains basic features of SOL plasmas, such as collisions, cross-field turbulent transport, parallel streaming, and particle loss at the parallel boundary, but without the magnetic gradient and curvature effects in a toroidal geometry. The configuration and parameters of the test case are similar to those in the experiments conducted by Carter and Maggs.²⁷ In the experiments, the LAPD creates a linear column plasma with a length L_z ∼ 17 m and a radius r ∼ 50 cm in the main chamber. The plasma is generated via ionization of neutrals by primary high-energy electrons injected by an anode–cathode arrangement from one end and is open on both ends along the magnetic field lines. The chamber wall can be biased with respect to the column plasma. For this study, we are concerned with the unbiased case; namely the wall is grounded.

There have been previous studies using various Braginskii-based fluid equations to study LAPD-type cases, by Rogers and Ricci,²⁸ Popovich et al.,^29,30 Friedman et al.,^31,32 and Fisher et al.³³ These looked at aspects of the turbulence characteristics and different possible driving mechanisms (including Kelvin-Helmholtz instabilities and nonlinear drift-wave mechanisms involving non-normal modes). Our simulation setup follows the fluid simulation by Rogers and Ricci,²⁸ and the subsequent gyrokinetic study by Shi et al.¹² These simulations produced results that have qualitative similarities to the experiments, though more detailed comparisons could be considered in future work, as discussed in Sec. V.

The remainder of this paper is organized as follows. Section II describes the basic equations. The parallel nonlinear term is included in the gyrokinetic equation for the gyrocenter distribution, and the Poisson equation is modified to the form used in Ref. 12. Numerical discretization schemes and the corresponding implementation of sheath boundary conditions are described in Sec. III. The setup and results of the simulation of an LAPD plasma are presented in Sec. IV. Section V contains conclusions and a discussion on future work.

By taking the electrostatic limit and the zero-Larmor-radius limit, the gyrocenter distribution in a straight field line geometry is evolved according to^22,23

where the E × B velocity is $vE=cB0b̂0×∇ϕ$⁠. In the global version of GENE, the gyrocenter distribution for species s is split into a static background depending on the x-coordinate and a fluctuation part $Fs(x,v∥,μ,t)=F0s(x,v∥,μ)+fs(x,v∥,μ,t)$⁠. Note that the last term is the parallel nonlinear term that is typically neglected for core plasmas due to the ordering $fs/F0s≪1$⁠.

For the purpose of consistent discretizations, we combine the parallel nonlinear term with the parallel linear term, the E × B nonlinear term with the E × B driving term, to obtain

Hereafter, the three terms on the right-hand side are termed the free-streaming term, the parallel nonlinear term, and the E × B nonlinear term, respectively. The free-streaming term combined with the parallel nonlinear term describes the dynamics parallel to the magnetic field lines, while the E × B nonlinear term describes the perpendicular dynamics.

By substituting the E × B velocity, introducing the field-aligned-coordinates with $B0=C∇x×∇y$ and $J−1=B0·∇z/C$⁠, and adding a collision term, the equation becomes

In this work, the Lenard–Bernstein collision operator^12,34 acting on full-f is applied

with $u∥,s=∫d3vv∥Fs/ns, msvT,ss′2=∫d3vms(v−us′)2Fs/3ns$⁠, $ns=∫d3vFs$⁠. The same-species collisions, including electron–electron (⁠$ss′=ee$⁠) and ion–ion (⁠$ss′=ii$⁠) collisions, thermalize the distribution function to a local Maxwellian, contain pitch-angle scattering and velocity-space diffusion, and conserve density, momentum, and energy. Standard expressions for the collision frequencies are used, $νee=42πneλe43meTe3/2$⁠, $νii=4πniλe43miTi3/2$⁠, where λ is the Coulomb logarithm (Huba 2016, p. 37).⁴⁶ The electron–ion (⁠$ss′=ei$⁠) collision term contains the collisional drag and pitch-angle scattering of electrons on ions. The electron–ion collision frequency is taken as $νei=νee/1.96$⁠, in order to approximately account for the parallel conductivity coefficient in a plasma. The ion–electron (⁠$ss′=ie$⁠) collision term is much smaller than the ion–ion collision term (⁠$νie/νii∼me/mi$⁠) and is neglected, as was also done in the Gkeyll study.¹² This leads to a small violation of the momentum conservation in the model.

By neglecting the Debye shielding, the nonlinear gyrokinetic Poisson equation in a straight field line geometry is^35–37 

where the long-wavelength limit $(k⊥ρi≪1)$ has been taken. The gyro-averaging is neglected in calculating the gyrocenter density, which is defined as

The potential ϕ is a nonlinear function of density because of the polarization charge density term. As an initial step for code developments, the nonlinear term is linearized by replacing the ion gyrocenter density $ni(x)$ with the background ion gyrocenter density n_i0, which is taken to be a constant in space and time, as done by Shi et al.¹² Thus, the linearized Poisson equation solved in the code is

We plan to implement a fast solver for the nonlinear Poisson equation (5) in the future [see a discussion on this in Sec. V].

The gyrokinetic equation (3) is discretized based on the so-called “method of lines.” A fourth-order explicit Runge–Kutta method is used for the time integrator. The derivatives in phase space contain three parts—the parallel dynamics, the E × B nonlinearity, and the collisions. In the global version of GENE, a Dirichlet or Neumann boundary and a fourth-order finite-difference method are applied in the x-direction to accommodate radial profile variations seen in the tokamak core, whereas a periodic boundary and a spectral method are used in the y-direction. The E × B nonlinearity is discretized with a mixed spectral/finite-difference variant of the Arakawa method.²³ It is implemented pseudospectrally in the y-direction and the so-called 3/2-dealiasing rule is applied.³⁸ For the LAPD test case, in addition to the E × B nonlinearity, different k_y modes interact also via the parallel nonlinearity, the open boundary, and the collisions. However, discretizations of these terms are not spectral and no dealiasing measures are taken. Therefore, for consistency and simplicity, the discretization of the E × B nonlinearity is modified to the finite-difference version of the Arakawa method in both x- and y-directions. In Fourier representation, this is achieved simply by replacing the multiplication factor ik_y with $i[8 sin (kyΔy)−sin (2kyΔy)]/(6Δy)$ for the first-order derivative, where Δy is the grid size. The former corresponds to discretizing the derivative spectrally, whereas the latter corresponds to discretizing the derivative with a fourth-order central difference stencil. In addition, the spectral multiplication factor for the second-order derivative in the Laplace operator of the field equation, $(iky)2$⁠, is replaced with a finite-difference one in Fourier representation, $[−cos (2kyΔy)+16 cos (kyΔy)−15]/(6Δy2)$⁠. After these two modifications, the code adopts a fourth-order finite-difference method consistently in the perpendicular directions. Finally, a small perpendicular numerical diffusion term is applied to allow for robust code operation.²² 

Below discretizations of the parallel dynamics, the parallel boundary and the collisions are described. We note that they are expressed in real space in both x- and y-coordinates; the forward and backward Fourier transformations along y-coordinate are applied wherever necessary.

The discretization of the dynamics in the parallel phase space $(z,v∥)$ is similar to that used in our previous work based on a 1D1V code.²⁵ For the purpose of describing sheath boundary conditions, we briefly summarize it.

First, the parallel dynamics is reformulated in a conservative form

where the free-streaming flux and parallel nonlinear flux are defined, respectively, as

and

Second, the spatial and velocity derivatives are discretized with a finite-volume method. Integrating Eq. (8) over each cell $Ik,l=[zk−12,zk+12]×[v∥,l−12,v∥,l+12]$ and applying the Stokes' theorem, we obtain the discretized equation

where $k=1,…,Nz, l=1,…,Nv$⁠. The cell-averaged distributions (indicated in the following with an overbar) are defined as

The free-streaming flux and the parallel nonlinear flux at the cell edges are defined, respectively, as

and

The net fluxes across the edges of cell I_k,l determine the change in $f¯s,k,l$⁠. By using the midpoint rule for the integration, the free-streaming flux and the parallel nonlinear flux are approximated as

and

Under this approximation, the overall order of the finite-volume method at most is $O(max(Δv∥,Δz)2)$⁠.³⁹ 

Third, the fluxes are reconstructed by first interpolating the data biased toward upwind direction with a piecewise polynomial, and then evaluating the polynomial at the cell edges. The free-streaming flux for $v∥>0$ at the edge $z=zk+12$ is

This results in a third-order upwind scheme for the free-streaming term

The species index s and the parallel velocity index l are omitted for clarity in Eqs. (17) and (18). The parallel nonlinear flux is similarly reconstructed. Its upwind direction is determined by the parallel derivative of the potential $∂ϕ/∂z$⁠, which is approximated by using a fourth-order central difference stencil for the inner points (⁠$3≤k≤Nz−2$⁠), a third-order difference stencil biased towards the center for the points next to the edge points (k = 2, N_z – 1), and a second-order difference stencil biased towards the center for the edge points (k = 1, N_z). These stencils are constructed with a standard technique,⁴⁰ namely, first interpolating the electrostatic potential inside the domain with piecewise polynomials, and then evaluating the derivatives of the polynomials at desired coordinates.

For typical edge and LAPD plasmas, the electron thermal velocity is much larger than the ion thermal velocity. In response to faster streaming to the wall by electrons than ions, a Debye sheath is formed at the plasma–wall interface to set up a potential drop from the plasma to the wall. The potential drop accelerates outgoing ions to the wall, slows down outgoing electrons, and reflects the low-energy electrons back to the plasma, so that the loss of electrons to the wall is approximately in balance with the loss of ions. However, the usage of the gyrokinetic Poisson equation allows the simulation spatial and temporal scales to be much larger than the sheath characteristic spatial and temporal scales, which are on the order of a few Debye lengths and plasma periods, namely $ρi≫λD, Ωi≪ωpe$⁠.

To account for the sheath effects on the plasma, a conducting-sheath boundary model is introduced by Shi et al.¹² In this model, the sheath potential is determined from the turbulence in the plasma, and the electron cutoff velocity v_c is calculated from the potential drop according to

where the sheath potential $ϕsh=ϕ(x,y,z=±π)$ at the upper $(zU=π)$ and lower $(zD=−π)$ boundaries is solved from the Poisson equation, and the wall is grounded $ϕwall=0$⁠. Outgoing ions are accelerated to the wall without reflection, while outgoing electrons with velocity lower than the cutoff velocity are repelled and reflected. Particle reflection at the upper boundary can be written as

where $zU+$ is the coordinate immediately outside the domain. The lower boundary conditions are similar. The cutoff velocity according to Eq. (19) typically falls within a cell $vc∈(v∥,l−12,v∥,l+12)$ and not exactly on a cell edge. The fraction of electrons reflected in the cutoff cell is approximated as

Physically, particles stream along the field lines and are lost to and reflected from the wall at the parallel boundary. In the finite-volume scheme, the sheath boundary enters the equation by setting the free-streaming flux for the edge cells. Define the last cell as $IN=[zN−12,zN+12]$⁠, then the domain edge (wall) is located at $zN+12$⁠. The upwind flux representing particles leaving from the system to the wall at the upper boundary $zU=zN+12$ is evaluated as

The second identity is valid because the first-order upwind flux is used for the last cell. The upwind flux representing particles reflected from the wall at $zU=zN+12$ is

where $Fs(zU+,v∥)$ for $v∥<0$ is determined from the corresponding outgoing flux (23) and the sheath reflection (20) and (21). In the code, $Fs(zU+,v∥)$ is stored in a boundary cell and updated in each stage of each Runge–Kutta time step.

In the conducting-sheath boundary model, net currents are allowed to flow through the wall, which is in contrast to the insulating logical-sheath boundary model used in previous studies.^25,41,42 In the logical-sheath boundary, the electron cutoff velocity is determined by requiring the loss of electrons to the wall to match the loss of ions, and the sheath potential is estimated from the cutoff velocity. Compared with the logical-sheath boundary, the conducting-sheath boundary reflects more properly the physical setup of the aforementioned LAPD experiments, in which the currents are free to flow along the field lines at both ends of the main chamber to adjust to the plasma state evolution self-consistently. In addition, the conducting-sheath boundary can handle the case of a biased wall by setting the wall potential ϕ_wall in Eq. (19).

The structure of the Lenard–Bernstein collision operator on full-f is similar to the test-particle part of the linearized Fokker–Planck collision operator on δf implemented in GENE. Therefore, we apply the same finite-volume method to discretize the outer derivative in the Lenard–Bernstein operator (4) to obtain

where the weights of the grids are defined as $Δμm=μm+12−μm−12$ and $Δv∥,l=v∥,l+12−v∥,l−12$⁠. A second-order central difference method is used for the inner derivative in Eq. (4) to evaluate the finite-volume fluxes across the edges of the cell $Il,m=[v∥,l−12,v∥,l+12]×[μm−12,μm+12]$ in velocity space

where $l=1,2,…,Nv, m=1,…,Nμ$⁠, and $αss′=vT,ss′2$ is introduced.

Similar to the implementation of the Fokker–Planck collision operator on δf, the implementation of the Lenard–Bernstein collision operator on full-f conserves the density by setting the fluxes at the boundaries. Specifically, the density change due to each type of collision can be expressed as

By setting $Ψv∥,12,m=Ψv∥,Nv+12,m=0$ for $m=1,2,…,Nμ$⁠, and $Ψv∥,l,12=Ψv∥,l,Nμ+12=0$ for $l=1,2,…,Nv$⁠, the density change is identically zero due to cancellations in the telescope sums.

The approach to conserving energy numerically is to treat $αss′$ as a free parameter, rather than directly evaluate it from its definition $αss′=vT,ss′2$ and use it in the finite-volume fluxes. This approach was introduced in previous gyrokinetic study in Ref. 12. The energy change due to each type of collision is

where the density conservation is used. By requiring the energy change to vanish, the parameter $αss′$ is solved for each type of collision.

We tested the implementation of the collision operator with a drifted-Maxwellian distribution and confirmed that the density and energy are conserved to machine precision, and that $αss′$ solved from the energy conservation condition is close to the value evaluated from its definition.

To test the extensions, we select the setup with LAPD-like parameters used in the gyrokinetic study by Shi et al.¹² The authors modified and adopted the setup from the earlier Braginskii-fluid study by Rogers and Ricci.²⁸ The nominal parameters are $Te0=6 eV, Ti0=1 eV, mi=3.973 mp$ (m_p is the proton mass), $B0=0.0398 T, n0=2×1018 m−3$⁠. The simulation domain is $(x,y,z)∈[−Lx/2,Lx/2]×[−Ly/2,Ly/2]×[−Lz/2,Lz/2]$⁠, with $Lz=18 m, Lx=Ly=L⊥=100ρs0=1.25 m$⁠, $ρs0=cs0/Ωi$⁠, and $cs0=Te0/mi$⁠.

For the LAPD-like parameters, the time steps allowed in the explicit Runge–Kutta integrator are determined by the electron–electron collisions. They are estimated from the diffusion term, namely $4νeevth2Δt/(Δv∥)2≤CCourant$⁠, where the Courant number C_Courant is taken to be 0.1. The time steps estimated from the collisions are typically smaller by a factor of ∼1000 than those constrained by the nonlinear terms. Furthermore, compared with previous gyrokinetic study by Shi et al.,¹² in which relatively low resolutions are allowed because the gyrokinetic equation is discretized with an energy-conserving discontinuous Galerkin method and the charge density integral is evaluated with the Legendre–Gauss quadrature, the GENE code requires higher perpendicular resolutions to resolve the perpendicular structures, and higher velocity-space resolutions to evaluate the collisions and the charge density. For a run with typical resolutions shown later, the small time steps combined with the relatively high resolutions would demand enormous computing resources if the experimental parameters were used. Two reductions are introduced to alleviate the problem. First, as done in Refs. 28 and 12, the electron-to-ion mass ratio is reduced to m_e/m_i = 1/400, so that the time steps are increased, while a sufficient separation of electron scales from ion scales remains. Second, in the gyrokinetic study by Shi et al.,¹² the electron–electron and electron–ion collision frequencies are reduced by a factor of 10 to increase the time steps. After the reduction, the electron mean-free-path is $Lmfp=vTe/νee≈0.57 m$ (using $ne=n0$ and $Te=3 eV$⁠); hence $Lmfp/Lz≪1$⁠. In our simulation, the electron–electron and electron–ion collision frequencies are reduced by a factor of 100. Correspondingly, the ion–ion collision frequency is reduced by a factor of 10. Hence, the electron mean-free-path is smaller than but comparable with the parallel length of the box; our simulation is marginally collisional.

To describe the initial conditions and the source term, we define $A(r;c)$ as a function that falls from the peak value of 1 at r = 0 to a constant value c for $r>L⊥/2$⁠,

where $r=x2+y2$⁠. The initial density profile for both ions and electrons is chosen to be $n0A(r;c=1/20)$⁠. The initial electron temperature profile has the form $5.7A(r;c=1/5)eV$ and the initial ion temperature profile is an uniform 1 eV. A top-hat-like source uniform along z representing the ionization of neutrals by primary electrons is applied continuously in time

where the radius of the source region is $rs=20ρs0=0.25 m,$ the spatial gradient scale length at the edge of the source region is $Ls=0.5ρs0=0.625 cm$⁠, and $FM,s(v∥,μ;Ts)$ represents a non-drift Maxwellian with a temperature T_s. The ion source temperature is $1 eV$⁠, and the electron source temperature is $6.8A(r;c=1/2.5) eV$⁠. The radial profiles of the electron source rate and temperature are shown in Fig. 1. The spatially integrated source rate is 6.5 × 10²¹ s⁻¹. This number is the total source rate used in the Gkeyll study,¹³ which includes a particle source of approximately 5.5 × 10²¹ s⁻¹ and a time-varying particle source of approximately 10²¹ s⁻¹ that comes from adding electrons and ions to the system in order to keep the density above a floor of 2 × 10¹² m⁻³ everywhere. The density is also kept above that floor level in GENE, but the rate of added source is negligible (below 0.1 × 10²¹ s⁻¹) throughout the simulation. Furthermore, it is checked and confirmed that the particle loss rate at the parallel boundary stays approximately at 6.5 × 10²¹ s⁻¹ after the system reaches a quasi-steady state.

In the global version of GENE, the background profiles of density and temperature are uniform along the y-coordinate, whereas in the present test case, both the initial profiles and the source term depend on the y-coordinate. Therefore, we set the background distributions to zero F_0s = 0 and f_s = F_s. The code calculates the full-f evolution.

We note that the numerical schemes in GENE are not designed to preserve the positivity of the distribution function. For the LAPD test case, we encountered the problem of negative distribution functions, which in principle can result from every term in the discretized version of the gyrokinetic equation. This problem was investigated in the Gkeyll study¹² (see Sec. 2.1.1 therein). The authors introduced a generic procedure to enforce the positivity of the distribution function. The cause and severity of the problem are not the same in the two codes, since the gyrokinetic equation is discretized with a discontinuous Galerkin method in Gkeyll and a finite-volume method in GENE. However, the positivity-adjustment procedure used in Gkeyll is also suitable for and implemented in GENE. The procedure consists of three steps and is summarized here. First, the negative values of the distribution functions are replaced with zeros. Second, to compensate for the increased density, the distribution functions are scaled down uniformly in velocity space to keep the density unchanged at fixed x for each species. Third, to remove the unphysical parallel energy introduced in the first two steps, a drag term in the parallel-velocity dimension is added to the gyrokinetic equation. The magnitude of the velocity drag is determined by requiring the final parallel energy at fixed $(x,μ)$ to match its value before the adjustment. Similarly, a drag term in the dimension of magnetic moment is implemented to correct the perpendicular energy. While the results to be shown in Sec. IV B appear physically reasonable, it is highly desirable to remove this positivity-adjustment procedure by discretizing the gyrokinetic equation systematically with a 5D finite-volume method and proper flux limiters so that the positivity of the distribution function is automatically preserved.⁴³ This requires substantial efforts and is subject to our future consideration.

Below we present the results from a typical run with the numbers of grid points N_x = 72, N_z = 32, N_v = 32, $Nμ=12$⁠, and y-coordinate wavenumber $Nky=36$ (equivalent to N_y = 72). The velocity-space domain is $(v∥,μ)∈[−4,4]2Ts,ref/ms×[0,16]Ts,ref/B0$⁠, with $Te,ref=3 eV$ and $Ti,ref=1 eV$⁠. The $v∥$-grid spacing is uniform, and μ-grid spacing is determined according to the Gauss–Laguerre quadrature so that the density integral achieves a sufficiently high accuracy. The minimum parallel and perpendicular energies represented by the velocity-space grids (assuming the distribution function collapses to the cell $I=[0,Δv∥]×[0,Δμ1]$⁠) are $T∥,e,min=0.05 eV$ and $T⊥,e,min=0.13 eV$ for electrons, and $T∥,i,min=0.017 eV$ and $T⊥,i,min=0.043 eV$ for ions. The run was performed on ∼4000 CPU cores for 24 h of wall-clock time. The total physical time of the run is ∼6 ms. For a reference, the ion transit time is $τi∼Lz/2/cs∼1.1 ms$⁠, using $Te=3 eV$ in calculating the ion acoustic velocity c_s. For this run, the background ion gyrocenter density in the linearized Poisson equation (7) was taken to be $ni0=6×1017 m−3$⁠.

Before presenting the results, we caution the readers that we do not claim our simulation which is a realistic representation of the experiments, for two major reasons. First, there are elements of the LAPD experiments that are not directly or fully modeled, such as the ionization and plasma heating process by the high-energy primary electrons, the ion–neutral collisions, and the parallel boundary conditions. Second, our gyrokinetic simulation is performed with a reduced mass ratio and a much lower collisionality to reduce the necessary computing resources. Lowering the collisionality from a highly collisional LAPD-like case to a marginally collisional one inevitably modifies the physical process. It is important to bear this in mind when comparing our simulation with previous fluid^28,33 and gyrokinetic¹² simulations. Our goal here is to show that the extended code produces physically reasonable features that are seen in experimental data and previous modeling results. Even though to a large extent our results agree with the Gkeyll results,¹² the focus of this paper is not benchmarking the two codes, which involves running a set of test cases with the same parameters and resolutions, preferably of low computation intensity, and quantifying and investigating the differences in simulation results as the parameters and resolutions are varied. This benchmark is important and necessary, but requires substantially more efforts, and will be conducted in the future.

Figure 2 shows the electron density evolution on the mid-plane. At the beginning of the simulation, the density builds up in response to the top-hat-like source, resulting in a steep density gradient at the edge of the source region. At $t≈700 μs$⁠, the plasma becomes unstable and the waves with an azimuthal integer wavenumber m ≈ 19 start to grow at the edge of the source region. As a consequence, a plasma rotation in the ion-diamagnetic-drift direction becomes visible in the density profile, as seen in previous fluid^28,33 and gyrokinetic simulations.¹² The waves are saturated progressively and turbulence develops and spreads radially, in both inward and outward directions. The system gradually reaches a quasi-steady state by $t≈3 ms$⁠.

To better visualize the wave generation and turbulence development, a movie showing the evolution of the electron density, electron temperature, and potential on the mid-plane is presented in the supplementary material. Figure 3 shows snapshots of the density, density fluctuations, temperature, and potential on the mid-plane after the turbulence reaches a quasi-steady state. The azimuthal structures near the edge of the source region in the temperature and potential profiles collocate with those in the density profile. All three quantities follow a similar trend of decreasing from the center outwards, although a small peak occurs in the potential near the edge of the source region. It is seen that little variation appears in the temperature and potential profiles at large radius, $r≳0.5 m$⁠, indicating that this low-density region is not significantly influenced by the turbulence and the perpendicular-boundary effects on turbulence are negligible. The density fluctuations will be discussed later. Figure 4 shows the same quantities at the same time as in Fig. 3, but the cut is made in the x – z plane at y = 0. The density, temperature, and potential peak near the center and decrease gradually towards the open boundaries. The filaments resulting from slicing the azimuthal structures seen in Fig. 3 extend almost uniformly along the field lines, indicating that the turbulence is highly anisotropic.

The time-averaged profiles as a function of the radius are presented in Fig. 5. The data are selected in the region of –4 m < z < 4 m, as done in Ref. 12. The considerations of selecting this region are that it is similar to the location where the measurements were taken in the LAPD experiments, and that there is little variation in the simulation along the field lines in this region. The time averaging is performed from t ≈ 3 ms to t ≈ 6 ms, with one sampling point per ∼1 μs (50 time steps). The time steps are nonuniform, and the gaps between sampling points are uneven. The time averaging for a quantity y is calculated by using a rectangular rule for the time integration, namely $y¯=∑i=0i=N−1yi(ti+1−ti)/(tN−t0)$⁠. The data are binned by radius. The features seen in the snapshots shown in Fig. 3 are confirmed in the averaged profiles. The density and temperature are almost constant at small radius and drop off significantly at middle radius. Compared with the top-hat-like source, the electron density profile is broadened by the turbulence at the edge of the source region. The potential profile shows a similar trend, and a small peak occurs near the edge of the source region. The turbulence influence appears to be confined in the region of $r≲0.5 m$⁠. Beyond this region, the density continues to drop off, whereas the temperature and potential reverse to larger values near the side wall. This is understandable considering that the source temperature near the side wall is 2.72 eV. The potential is smaller and the turbulence fluctuations are weaker when compared with the Gkeyll results.¹² The temperature and potential drop off monotonically from the center to the side wall in that study. Note that the uniform μ-grid spacing restricts the temperature in that particular simulation to above 1.2 eV. The fluid simulation results^28,33 show that the temperature falls to close to zero near the side wall, but the source temperature in their simulations vanishes at large radius.

The statistics of density fluctuations are presented in Fig. 6. The density fluctuation is defined as $ñe(x,t)=ne(x,t)−n¯e(x)$⁠, where the time-averaged density is subtracted. The data and the time averaging process are the same as those used in Fig. 5. The root-mean-square (RMS) of the density fluctuations is calculated and normalized by the peak time-averaged density at r = 0. The strongest fluctuations occur at the edge of the source region, and the maximum fluctuation level is ∼13%. The amplitude of the fluctuations decays both radially inward and outward from the peak location. However, when normalized to the local background density, the peak fluctuations are at the 40% level and occur outside of the source edge. The power spectral density of the fluctuations is calculated by using a Fourier transformation and normalized to the value at 1 kHz. The data are the same as those in Fig. 5, but are further restricted to the region of 0.25 m < r < 0.3 m, and are interpolated linearly in time before the Fourier transformation is applied. The power spectral density shows that the fluctuations are broadband below ∼20 kHz and obey a powerlaw distribution with a powerlaw index n = –1.87 between 20 kHz and 400 kHz. Here the upper limit of this frequency range is set by the sampling frequency. The density fluctuation statistics agree with the Gkeyll results,¹² suggesting that the turbulence nature is the same in the two gyrokinetic studies. The previous fluid study³³ argues that the density fluctuations follow an exponential distribution in the frequency range between 4 kHz and 20 kHz.

The probability density function (PDF) of the density fluctuations has been used as an indicator of density blobs (structures of density enhancement) and holes (structures of density depletion) in LAPD turbulence.⁴⁴ A negatively skewed PDF indicates presence of density holes, whereas a positively skewed PDF indicates density blobs resulting from the plasma radial transport from the high-density source region to the low-density region in LAPD. The PDFs at three radial locations are calculated and presented in Fig. 7. The PDF is negatively skewed inside of the edge of the top-hat-like source region. It is close to a symmetric Gaussian distribution near the edge of the source region and turns positively skewed outside of the source region. While this trend is consistent with experimental data,⁴⁴ and previous fluid³³ and gyrokinetic¹² results, large-amplitude fluctuation events near the edge of the source region occur less frequent than a Gaussian distribution. A snapshot of the density blobs and holes is shown in Fig. 3. They appear side-by-side with each other and are results of the turbulent fluctuations.

In addition to the diagnostics that are commonly used for comparisons with experimental data, we have also examined the distribution functions, to further check the implementation of the open boundary and collisions, and to understand the effects of particle loss at the open boundary. Figure 8 shows the diagnostics data at the upper open boundary taken at the same physical time as the mid-plane snapshots shown in Fig. 3. Similar to the profiles on the mid-plane, the azimuthal structures in the density profile approximately collocate with those in the potential profile. As expected, the peak density and potential are smaller than those on the mid-plane. The electron distribution function shown in Fig. 8(c) is sampled near the center of the perpendicular plane and is representative of distribution function in the high-density region. The distribution function is normalized by $ns,ref/vTs,ref3$⁠, where the reference density is $ns,ref=2×1018 m−3$ for both species, and the reference thermal velocity is $vTs=2Ts,ref/ms$ with $Te,ref=3 eV$ and $Ti,ref=1 eV$⁠. The $v∥$-coordinate of each species is normalized by its reference thermal velocity, and the μ-coordinate is normalized by $Ts,ref/B0$⁠. For a reference, the drifted-Maxwellian distribution [Fig. 8(d)] is shown along with the data. The drifted-Maxwellian is constructed by using velocity moments (density, parallel velocity, and temperature) of the data. By visual inspection, the electron distribution is close to the drifted-Maxwellian. In the sheath boundary model, outgoing electrons with sufficiently high parallel-energies can overcome the electrostatic potential and stream to the wall, whereas lower-energy electrons are repelled and reflected back to the system, with their velocities reversed. This selective loss and reflection is expected to cause velocity space depletion. However, the diffusion term in the collisions can shuffle electrons in velocity space to fill in the depletion. When the distribution function is sliced and compared with the drifted-Maxwellian at fixed μ as shown in Fig. 8(e), the depletion effect at velocities greater than the cutoff velocity becomes clearly visible. By checking the distribution functions at different locations along the field line, we found that the depletion effect becomes weaker and gradually disappears as the probe is moved farther away from the boundary. For example, the distribution becomes nearly symmetric about $v∥=0$ at z = 6 m [Fig. 8(g)]. The positive shift is not visible in the normalized electron distribution because the electron reference thermal velocity is much larger than the drift velocity, but it is clearly seen in the ion distribution function at the boundary [Fig. 8(i)]. The depletion effect is not present in the ion distribution function probably because the streaming loss of ions is too weak to compete with the collision diffusion. Finally, we note that the distribution functions appear to be well-resolved in $v∥$⁠, and improvements are expected by increasing the resolution in μ to better resolve the high-μ portion, as indicated by Figs. 8(f), 8(h), and 8(j).

We have reported our recent progress regarding the extensions of the global version of the 3D2V gyrokinetic code GENE for applications of electrostatic edge plasma turbulence. The code employed a δf-splitting technique and was originally designed for simulating core turbulence. To accommodate the large-amplitude fluctuations that are often seen in edge plasmas, the parallel nonlinear term is included and discretized along with the corresponding linear term by using a finite-volume method. By including the parallel nonlinearity, the δf code becomes equivalent to and can be run as a full-f code. To account for sheath effects that are important in open-field-line systems such as the SOL, a conducting-sheath boundary model is implemented in the context of the finite-volume method. The sheath boundary model takes into account particle loss to and reflection from the wall, without directly resolving the sheath details. Finally, the Lenard–Bernstein collision operator acting on full-f is implemented. The implementation conserves the number density and energy. The extended code can simulate collisional electrostatic turbulence in straight open-field-line systems.

To test these extensions, the code is applied to simulating a test case based on LAPD experiments.²⁷ The test case has been modeled previously by Braginskii-fluid^28,33 and gyrokinetic¹² codes. Our simulation produces features that are seen in the experimental data and previous simulations. Turbulence develops and spreads radially inward and outward after the waves grow at the maximum density gradient region and are saturated progressively in time. The absolute amplitude of density fluctuations peaks at the edge of the top-hat-like source region, and there the density profile is flattened by the turbulence. The density fluctuations of the fully developed turbulence are broadband below ∼20 kHz and follow a powerlaw distribution in the higher-frequency range. Density holes are more likely present inside of the source edge, while density blobs are more likely present outside. The depletion of the electron distribution functions caused by the electron loss to the wall is largely filled due to collision diffusion in velocity space, and the distribution functions are close to the drifted-Maxwellian near the open boundary. For more detailed comparisons with LAPD experiments, future extensions could include effects that may alter the predicted rotation profiles, and thus the turbulence characteristics. These include the effects of the mesh anode, energetic electron emission from the hot cathode, and viscosity. In addition, it is highly desirable to implement the more sophisticated Fokker-Planck operator to conserve particle, momentum, and energy numerically,⁴⁵ preferably with an implicit time-stepping scheme, so that the collision effects on LAPD and edge plasmas can be realistically evaluated. Quantitatively matching GENE simulations with experimental results and examining the transport mechanisms in LAPD turbulence^28–33 are of great interest and left for future work.

A few additional improvements and complexities can be built up upon the present code. First, the Poisson equation (7) uses a constant background ion density in the polarization charge density term. To evaluate the effects of this simplification, a simulation with the background ion density doubled to n_i0 = 1.2 × 10¹⁸ m⁻³ has also been performed. No material change has been found in the averaged potential profile. This is understandable because the averaged potential at the sheath entrance needs to satisfy the condition ϕ_sh ∼ 3T_e to reflect most of the electrons to give an electron flux close to the ion flux to the wall; changing n_i0 will just change the small difference between the electron density and the ion guiding center density that is needed to maintain the sheath potential. In addition, it is found that the maximum fluctuation level in the potential is reduced slightly from ∼16% to ∼14%. Although the turbulence appears to be insensitive to the change in n_i0, we plan to develop a field solver to use a spatially dependent ion density in a quasi-steady state, or update the ion density in every time step to remove this limitation. The latter approach requires a fast iterative solver based on the conjugate gradient method to solve the resultant sparse-matrix equation from the nonlinear Poisson equation. Second, to simulate SOL turbulence, it is necessary to include the effects of magnetic gradient and curvature (see Chap. 5 in Ref. 13). This requires implementation of additional terms in the parallel nonlinearity that are neglected. Similar to the E × B nonlinearity and the parallel nonlinearity in the straight-field-line geometry, the nonlinear terms associated with the magnetic gradient and curvature can be combined with the linear terms and discretized by using a finite-volume method. Third, physical effects such as electromagnetic effects, X-point geometry, and plasma–wall interactions are required to realistically model the edge plasmas. A few of these directions are being pursued, and the progress will be reported in the future.

See supplementary material for a movie showing the evolution of the electron density n_e (in 10¹⁸ m⁻³), electron temperature T_e (in eV), and electrostatic potential ϕ (in V) on the mid-plane.

We thank the referee for suggesting examining the distribution functions at the parallel boundary. We thank T. Carter for useful discussions about LAPD, and B. Rogers and P. Ricci for discussing the test problem. Q. Pan would like to acknowledge G. Merlo and M. Weidl for helpful discussions. E. Shi would like to acknowledge T. Stoltzfus-Dueck for helpful discussions. This work was funded in part by U.S. Department of Energy under Contract No. DE-AC02-09CH11466 at the Princeton Plasma Physics Laboratory, through the Center for the Study of Plasma Microturbulence SciDAC project and the Max-Planck-Princeton Center for Plasma Physics. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Gkeyll simulations in support of this work were performed on the Perseus cluster at the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology's Research Computing department. Part of this work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under Grant Agreement No. 633053. Simulations have been performed on the EUROfusion High Performance Computer (Marconi-Fusion). The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s704.