Full-f gyrokinetic simulation of turbulence in a helical open-field-line plasma

Satisfactory and reliable quantitative predictions of turbulence and transport in the tokamak edge and scrape-off-layer (SOL) regions are widely believed to require the use of expensive gyrokinetic simulations in some capacity.^1–4 Some major outstanding questions that require numerical investigation include how the SOL power width is set,^5–7 how a confined plasma transitions from a low-confinement L mode to a high-confinement H mode,^8,9 and how high the H-mode pedestal temperature can get since the pedestal temperature has a major impact on the core temperature profile and the resulting gain.^10,11 Gyrokinetic simulations in the edge and SOL regions are challenging for several reasons (e.g., large-amplitude fluctuations, steep profile gradients, closed and open magnetic field lines, X-point, and effective sheath-model boundary conditions), but specialized particle-in-cell^12,13 and continuum gyrokinetic codes^14–18 have been making steady progress towards the ultimate goal as a predictive tool for boundary-plasma modeling. We refer the reader to Cohen and Xu¹ for a summary of early approaches to the gyrokinetic edge and SOL simulation and Krommes¹⁹ for an introduction to gyrokinetics. The particle-in-cell-based XGC1 code²⁰ is the only gyrokinetic code at present that is able to simulate turbulence in a three-dimensional diverted geometry. The first gyrokinetic simulations using continuum algorithms to simulate turbulence on straight open field lines were presented in Shi et al.¹⁴ and then in Pan et al.¹⁶ Here and in Shi,¹⁵ we extend this earlier work to present the first gyrokinetic continuum simulations on open field lines including curved toroidal fields, which can strongly enhance the drive of plasma instabilities.

The SOL refers to the tokamak plasma region of open magnetic field lines between the last closed flux surface (LCFS) and the first wall. Here, the field lines intersect material surfaces that act as plasma sinks where the loss rate of electrons and ions is kept in approximate particle balance by a Debye sheath layer. Plasma–surface interactions²¹ at the material interfaces can also contaminate the plasma with wall materials, which can severely degrade the fusion-plasma quality, but we do not yet incorporate these effects in the model being presented here.

Probe and imaging diagnostics have revealed the existence of intermittent coherent structures in the SOL referred to as plasma filaments or blobs,^22–26 which convectively transport particles, heat, momentum, and current across magnetic field lines.²⁷ Blobs are characterized by densities that are much higher than local background levels, a structure that is highly elongated along the magnetic field (much larger than the plasma minor radius) and much smaller scales perpendicular to the magnetic field, ∼10ρ_i, where ρ_i is the ion gyroradius.^27,28 Cross-field transport in the far SOL is highly intermittent due to blob propagation²⁸ and is consequently poorly described in terms of effective diffusion coefficients and convective velocities.²⁹ 

In a tokamak, the curvature and $∇B$ forces are believed to set up a charge-separated dipole potential structure across the blob cross-section which results in its outward radial propagation via convective $E×B$ transport.^27,30 Finite-temperature effects of the blob can also cause spin motion if the blob is sheath-connected, which can reduce this radial motion.³¹ Numerically, blobs' dynamics have been studied using seeded-blob fluid simulations.^32–35 Self-consistent blob formation has been studied with two-dimensional models^36–38 and in three-dimensional turbulence simulations.^12,39–41

The work presented here builds on our previous efforts in simulating open-magnetic-field-line turbulence in the Large Plasma Device⁴² using the gyrokinetic continuum capabilities of the Gkeyll code.¹⁴ In that study, the magnetic field was straight and uniform, and the plasma was highly collisional, which necessitated the use of an artificial electron-to-ion mass ratio (m_i/m_e = 400) and reduced electron collision frequencies to make the simulations tractable, given the explicit algorithm used at present for the collision operator. Nevertheless, we found that our numerical approach based on discontinuous Galerkin methods and sheath-model boundary conditions for an open-field-line region was stable and produced qualitatively reasonable results, which led to the first demonstration of open-field-line turbulence with a gyrokinetic continuum code. The reduced electron mass and collision-frequency restrictions have been relaxed for the simulations presented in this paper, which now also include a more sophisticated magnetic geometry.

We have added magnetic curvature and $∇B$ drifts to the Gkeyll code and can simulate a helical magnetic geometry approximating that in simple magnetized tori (SMTs), such as TORPEX⁴³ and Helimak.⁴⁴ In contrast to early work on the simulation of turbulence in SMTs based on the drift-reduced Braginskii equations and neglecting the ion temperature,^45–47 the gyrokinetic approach can naturally investigate plasmas with $Ti≳Te$⁠, which is commonly observed in the SOL.^48–50 Recent fluid simulations also include finite T_i.^51,52

Although our simulations do not yet simultaneously contain open- and closed-field-line regions,^51–56 we believe that many basic properties of SOL turbulence and transport are reproduced in this model. Additionally, the turbulence in this helical open-field-line geometry with parameters appropriate for a tokamak SOL has not been previously studied using a gyrokinetic PIC approach, either. We do acknowledge, however, that gyrokinetic PIC codes that have the necessary capabilities for the problem described in this paper have already been developed,^12,13 and it should be straightforward for these codes to implement this simple helical geometry for cross-code comparisons.

We discuss the details of the helical-SOL model in Sec. II, including equations solved, simulation geometry, and boundary conditions. Additional details about the underlying algorithms can be found in Refs. 14 and 15. We present simulation results obtained using the Gkeyll code in Sec. III, such as heat-flux profiles, fluctuation statistics, and particle fluxes. Our conclusions are given in Sec. IV. To facilitate future code comparisons, we also present the details of the initial conditions used in our simulations in the  Appendix.

In the non-orthogonal field-aligned geometry used in the simulation, z measures distances along field lines relative to the midplane (poloidal angle θ_pol = 0 in a tokamak), x is the radial coordinate, and y is constant along a field line and measures distances perpendicular to field lines. The simulation geometry is a flux tube on the outboard side that wraps around the torus a specified number of times, terminating on material surfaces at each end in z. The resulting mapping from field-aligned coordinates (x, y, and z) to standard cylindrical coordinates (R, φ, and Z) is given by R = x, $Z=z sin θ$⁠, and $φ=(y sin θ+z cos θ)/Rc$⁠, where R_c = R₀ + a and the field-line pitch $sin θ=Bv/B$ are taken to be constant, R₀ is the device major radius, a is the device minor radius, and B_v is the vertical (or poloidal) magnetic field. This simple helical geometry has vertical flux surfaces (i.e., ignores flux expansion) and has no magnetic shear, and some further approximations to differential operators are made assuming short-wavelength turbulence for now. The final model nevertheless includes the main effect of the bad-curvature drive by toroidal magnetic fields while using an efficient field-aligned grid. See Refs. 57–59 and 15 for further details.

We solve a full-f gyrokinetic equation written in the conservative form^60–62 

where $fs=fs(R,v∥,μ,t)$ is the gyrocenter distribution function for species s, $J=B∥*$ is the Jacobian of the gyrocenter coordinates, $B∥*=b·B*, B*=B+(Bv∥/Ωs)∇×b$⁠, and C[f_s] represent the effects of collisions, Ω_s = q_sB/m_s, and $Ss=Ss(R,v∥,μ,t)$ represents plasma sources. The characteristics are calculated as $Ṙ={R,H}$ and $v̇∥={v∥,H}$⁠, where the gyrokinetic Poisson bracket operator is

and the gyrocenter Hamiltonian is $Hs=12msv∥2+μB+qsϕ$⁠, where the long-wavelength limit has been taken to neglect gyroaveraging. A conservative Lenard–Bernstein collision operator⁶³ that neglects the velocity dependence of the collision frequency is used to model self-species and electron–ion collisions. Expressions for velocity-independent collision frequencies are taken from the NRL Plasma Formulary⁶⁴ [see Eqs. (2.52) and (2.53) of Ref. 15 for the exact expressions used].

This system is closed by the long-wavelength gyrokinetic Poisson equation with a linearized ion polarization density

where $ρs0=cs0/Ωi, cs0=Te0/mi$⁠, and $ni0g$ is the background ion gyrocenter density that we take to be a constant in space and in time. The long-wavelength approximations and the use of a linearized ion polarization density are made in a consistent manner, and so, the system described by Eqs. (1)–(3) conserves the number, momentum, and total energy. This system can be derived using the Lagrangian field theory⁶¹ by making all approximations at the system-Lagrangian level (e.g., see Appendix A of Ref. 65 for details).

In these equations, we neglect geometrical factors arising from a cylindrical coordinate system everywhere except in $B*=B+(Bv∥/Ωs)∇×b$⁠, where we make the approximation that perpendicular gradients are much stronger than parallel gradients

Here, we assume that $(∇×b)·ey=−1/x$⁠, where $ey=∇y$ is a “co-basis” direction. This type of approximation has also been employed in some fluid simulations of SMTs.^46,66 We assume that $B=Baxis(R0/x)ez$⁠.

Periodic boundary conditions are applied to both f and ϕ in y, and the Dirichlet boundary condition ϕ = 0 is applied in x, which prevents gyrocenters from crossing the surfaces in x. Conducting-sheath boundary conditions are applied to f in z, which partially reflect gyrocenters of one species and fully absorb gyrocenters of the other species into the wall depending on the sign of the sheath potential. The potential is determined by solving the gyrokinetic Poisson equation [Eq. (3)]. Evaluating this potential at the sheath entrances (the ends of the simulation domain in z) gives the sheath potential, which is used to determine which particles are reflected by the sheath. This is the gyrokinetic analog of how fluid codes have used the vorticity to calculate the potential and sheath effects (for example, see Ricci and Rogers⁴⁶). We refer to these boundary conditions as conducting-sheath boundary conditions^14,15 because they allow self-consistent currents locally in and out of the end plates. This is in contrast to the logical-sheath boundary-condition model,^67–69 which assumes an insulating sheath with zero current density at the end plates everywhere. There is no closed-field-line region in our present model.

We use parameters roughly approximating a singly ionized H-mode deuterium plasma in the NSTX SOL:^70,71 $ni0g=7× 1018 m−3$⁠, T_e ∼40 eV, T_i ∼60 eV, B_axis = 0.5 T, R₀ = 0.85 m, and a₀ = 0.5 m. Although we use parameters for an H-mode plasma, we do not attempt or claim to capture H-mode physics (e.g., an edge transport barrier) in our simulations since they include only the SOL and not the pedestal.

The simulation box has dimensions L_x = 50ρ_s0 ≈ 14.6 cm, L_y = 100ρ_s0 ≈ 29.1 cm, and $Lz=Lp/ sin θ$⁠, where L_p = 2.4 m, ρ_s0 = c_s0/Ω_i, and θ is the magnetic-field-line incidence angle. The magnetic field is taken to be composed primarily of a toroidal component with a smaller vertical component (referred to as B_v), resulting in a helical-field-line geometry that roughly approximates the tokamak SOL. We present results with $sin θ=Bv/Bz=(0.2,0.3,0.6)$ in Sec. III, which correspond to L_z = (12, 8, 4) m. The connection length to the divertor plate in the real NSTX experiment is typically quite long, over 10 m, but we consider smaller values that might represent the shorter connection length from the midplane to the X-point region, where the magnetic shear is very strong. In this study, the magnetic-field-line incidence angle is not accounted for in the sheath boundary conditions (i.e., no Chodura sheath⁷²).

We use an energy-conserving discontinuous Galerkin method for the spatial discretization of the equations, which is a generalization of the algorithm of Liu and Shu⁷³ for two-dimensional incompressible flow in the vorticity–stream function formulation. Time discretization is performed using an explicit third-order strong-stability-preserving Runge–Kutta algorithm.⁷⁴ The positivity of the distribution function is not automatically guaranteed in our algorithms, and our method to keep f > 0 results in the addition of a small amount of numerical heating ∼10% of the source power to the system. The details of the numerical algorithms, energy conservation, and sheath boundary conditions are discussed in Refs. 14 and 15.

The position-space extents are $x∈[R0+a0−Lx/2,R0+ a0+Lx/2], y∈[−Ly/2,Ly/2], z∈[−Lz/2,Lz/2]$⁠, and the velocity-space extents are $v∥s∈[−v∥s,max,v∥s,max]$⁠, where $v∥s,max=4vts=4Ts/ms$ and $μs,max=(3/4)msv∥s,max2/(2B0)$⁠, where B₀ = B_axis R₀/(R₀ + a₀). The solution in each cell is expanded using piecewise-linear basis functions, i.e., the span of monomials in the five phase-space variables with each variable degree ≤1. This choice results in 32 degrees-of-freedom per element to represent the distribution function and Hamiltonian. The grid resolution is $(Nx,Ny,Nz,Nv∥,Nμ)=(18,36,10,10,5)$⁠, and a uniform grid spacing is used.

The plasma density source has the following form:

where x_S = –0.05 m + R₀ + a₀, λ_S = 5 × 10⁻³ m, and S₀ is chosen so that the source has total (electron plus ion) power P_source = 0.27L_z/L_z0 MW, where L_z0 = 4 m. The expression for the source power comes from multiplying P_SOL = 5.4 MW, the total power into the SOL, by the fraction of the total device volume covered by the simulation box. A floor of 0.1S₀ is used in the $|z|<Lz/4$ region to prevent regions of n ≪ n₀ from developing at large x, which can result in distribution-function positivity issues. The distribution function of the sources is non-drifting Maxwellians with a temperature profile T_e,i = 74 eV for x < x_S + 3λ_S and T_e,i = 33 eV for x ≥ x_S + 3λ_S. These choices result in an integrated source particle rate of ≈9.6 × 10²¹ s⁻¹ for the L_z = L_z0 (B_v/B_z = 0.6) case.

We do not yet include a closed-field-line region in our simulations, so we only simulate a SOL. The x < x_S + 3λ_S region will be referred to as the source region in this paper, while the x ≥ x_S + 3λ_S region will be referred to as the SOL region. We can think of the x = x_S + 3λ_S location as the LCFS.

Starting from an initial condition estimated by the steady-state solution of one-dimensional fluid equations (see the  Appendix), the sources steepen the plasma profiles, quickly triggering curvature-driven modes that grow on a timescale comparable to $γ∼cs/Rλp$⁠. We emphasize that our system does not contain ballooning modes since there are no “good-curvature” regions. As shown in Fig. 1, radially elongated structures extending far from the source region are generated and subsequently broken up by sheared flows in the y direction in the source region, leaving radially propagating blobs in the SOL region. Using the time-averaged profiles from the same B_v/B_z = 0.3 (L_z = 8 m) simulation, we estimate γ ∼1.9 × 10⁵ s⁻¹ using λ_p ≈ 2.9 cm, T_e ≈ 30 eV, and R = x_S = 1.3 m. On a time scale long compared to γ⁻¹ and τ_i = (L_z/2)/v_ti ∼50 μs, the conducting-sheath boundary conditions maintain a quasi-steady state in which the particle losses to the end plates are balanced by the plasma sources. Snapshots of the electron density, electron temperature, and electrostatic potential from the quasi-steady state (t = 625 μs) for the B_v/B_z = 0.3 case are shown in Fig. 2.

For the steepest magnetic-field-line-pitch case (B_v/B_z = 0.6), we performed a second simulation with magnetic-curvature effects removed, keeping all other parameters unchanged. The resulting magnetic geometry consists only of straight magnetic field lines, so coherent structures of elevated plasma density cannot become polarized by curvature forces. As shown in the electron-density snapshot comparison in Fig. 3, the presence of magnetic curvature appears to have an important role in the turbulent dynamics of the system. When magnetic-curvature effects are removed, the radial propagation of coherent structures into the SOL region appears to be significantly weakened or absent, and most of the density is localized to the source region.

Figure 4 compares radial profiles of the background electron densities, normalized (to the local background value) electron-density fluctuation levels, and radial E × B particle fluxes Γ_n,r between these two simulations. The radial particle flux due to electrostatic turbulence is calculated as $Γn,r=⟨ñeṽr⟩$⁠,²⁸ where v_r = E_y/B and the brackets $⟨…⟩$ indicate an average over a period that is long compared to the fluctuation time scale and an average over y and the central region in z, –0.5 m < z < 0.5 m. The fluctuation of a time-varying quantity A is denoted as $Ã$⁠, which is related to the total A as $Ã=A−⟨A⟩t$⁠. Here, the brackets $⟨…⟩t$ indicate an average in time. Notable differences between these two simulations are found in all three quantities plotted. Compared to the helical-SOL simulation, the straight-field-line simulation has a background density profile that decays more rapidly, fluctuation levels that quickly drop to ≈0% outside x ≈ 1.35 m, and a ≈ 2.5 times smaller Γ_n,r that also drops to approximately zero outside x ≈ 1.34 m.

The 10%–30% relative fluctuation levels in n_e shown in Fig. 4 call into question the use of a linearized ion polarization density in the gyrokinetic Poisson equation [Eq. (3)]. Additionally, the spatial variation of the background n_e profile calls into question the use of a spatially constant $ni0g$ in the same term. The code currently does not have the option of using a nonlinear ion polarization density (which also requires a higher-order term in the Hamiltonian for consistency), so we leave an investigation of this approximation to future work. Although the results indicate that the approximations made to the ion polarization density in Eq. (3) are not formally valid, it is likely that the use of a nonlinear ion polarization density will (slightly) change σ_g instead of ϕ. As we show later, it appears that the potential in these simulations is primarily determined by the parallel force balance in the electrons.

We have also performed a scan of the mass ratio m_i/m_e from the actual ratio of 3698 down to 100 (by increasing the electron mass), and we observed no significant quantitative or qualitative changes in the turbulence. The mass ratio might play an important role in a different parameter regime, however.

Effects connected to B_v ∼ B_p ∼ I_plasma are explored by changing the magnetic-field-line incidence angle since $sin θ=Bv/Bz$⁠. We have performed simulations at three values of magnetic-field-line pitches B_v/B_z = (0.2, 0.3, and 0.6), which correspond to L_z = (12, 8, and 4) m and θ = (20.14∘, 30.47∘, and 64.4∘). We scale the source appropriately in each simulation to keep the volumetric source rate the same. In all these simulations, the source is localized to the z ∈ [–L_z/4, L_z/4] region to model a source with a fixed poloidal extent. As θ is decreased, the plasma profiles are observed to become less peaked, implying that turbulence transport in the x-direction increased with decreasing θ.

We calculate the steady-state parallel heat flux $q=∑s∫d3v Hsv∥fs$ at the sheath entrance and average q in the y-direction to obtain a radial profile of the steady-state parallel heat flux for each case. To compare the heat fluxes on an equal footing, we plot the component of the parallel heat flux normal to the divertor plate $q⊥=q sin θ$ in Fig. 5. Compared to the B_v/B_z = 0.6 case, the heat-flux profiles for the cases with a shallower pitch are much broader. This behavior is consistent with the observation in tokamaks that the SOL heat-flux width is inversely proportional to the poloidal magnetic field (analogous to B_v in this model) and the plasma current,^5,75 although the physical reasons behind the scaling in our model and in a tokamak SOL may be quite different. We note that a significant amount of plasma in the smallest θ case gets near the outer radial wall, where further radial transport is suppressed, since the outer boundary is taken to be an ideal conducting plate with constant ϕ, so the E × B velocity into the side walls, ∝ ∂ϕ/∂y, vanishes. Simulations with a larger domain extent in the x coordinate (and/or finite-Larmor-radius effects in the collision operator to include classical transport to the side wall) might exhibit more of an exponential fall off over a wider radial range, further reducing the density in the right-hand side of the simulation.

The broad heat-flux profiles in Fig. 5 can be connected to the increased outward radial turbulent transport as B_v/B_z becomes shallower. We compute the steady state radial particle flux Γ_n,r near the midplane in the region –0.5 m < z < 0.5 m for each value of B_v/B_zx and plot the y-averaged fluxes in Fig. 6 (solid lines). Since the simulation box occupies a larger fraction of the device volume as B_v/B_z is decreased, but the source occupies the same fraction of the simulation box and has a fixed volumetric source rate, the background density levels increase as B_v/B_z decreases. Another way to say this is that with a fixed volumetric source density (fixed in particles per cubic meter per second), the mean density is expected to increase as the parallel connection length L_z/2 increases and the parallel loss rate ∼1/τ_i = 2v_ti/L_z decreases. Therefore, the magnitude of the Γ_n,r profiles in Fig. 6 should not be taken alone as a measure of turbulence levels.

The Γ_n,r profiles can be compared with the radial particle fluxes that result from assuming Bohm diffusion, i.e., Γ_B = D_B∂_xn_e, where the diffusion coefficient D_B = (1/16)k_BT_e/(eB). In the x > 1.36 m region, Γ_n,r/Γ_B ≈ 16 for the B_v/B_z = 0.2 case, while Γ_n,r/Γ_B ≈ 8 for the B_v/B_z = 0.6 case. One might expect the maximum level of turbulent transport to be comparable to the levels set by D_B, but it is important to remember that D_B is a diffusive transport estimate. The convective transport of blobs in these simulations appears to be responsible for the much-higher turbulent fluxes. Experimental data from tokamaks also suggest that the higher-than-Bohm particle transport in the SOL is due to the non-diffusive transport of blobs.^28,76

Density fluctuation statistics are often of interest in the SOL to characterize the turbulence. Considering again a time-varying quantity A, we define the skewness of A as $E[Ã3]/σ3$ and the excess kurtosis of A as $E[Ã4]/σ4−3$⁠, where σ is the standard deviation of A and E[…] denotes the expected value. Figure 7 shows the radial profiles of the normalized fluctuation level, skewness, and excess kurtosis for electron-density fluctuations and electrostatic-potential fluctuations computed near the z = 0 m plane. The density and potential fluctuations are normalized to their local background values. The positive skewness and excess kurtosis values are signatures of intermittency, which indicates an enhancement of large-amplitude positive-density-fluctuation events and is connected to the transport of blobs.^28,77

A somewhat counter-intuitive result is the reduction of density fluctuation levels as B_v/B_z is decreased, given that Figs. 5 and 6 indicate that turbulent spreading is increased as B_v/B_z is decreased. The skewness and excess kurtosis plots in Fig. 7 indicate that the density fluctuations become closer to a normal distribution as B_v/B_z is decreased. These trends in the density fluctuation statistics can be understood by noting that the background density profile becomes less peaked and more uniform in the x-direction as B_v/B_z is decreased, so a blob that is formed in the source region propagating in the SOL has a density that is closer to the background level, which results in lower relative fluctuation, skewness, and excess kurtosis values when compared to the large B_v/B_z case. Additionally, the density flux is constrained by the use of a fixed volumetric source rate, so as the background density increases with decreasing B_v/B_z, the relative density fluctuation levels tend to decrease. We also observe that the potential fluctuations are much less intermittent than the density fluctuations at the same B_v/B_z. This observation could be a real, physical effect, but we note that the fact that the temperature at large x runs into the grid resolution (the lowest temperature that can be represented on the velocity grid) could be influencing the potential fluctuation statistics in this region. Unlike the density fluctuations, the normalized potential fluctuation levels tend to increase with the decreasing B_v/B_z.

Figure 8 shows the radial profiles of the steady-state ion and electron temperatures and ion-to-electron temperature ratios near the midplane for different B_v/B_z values. For all three simulations, T_i/T_e falls in the range of 1.5–2, which is within the range of 1–10 that is observed a few centimeters outside the LCFS in tokamaks.⁴⁹ Similar to the heat-flux profiles shown in Fig. 5, the profiles are the steepest for the case with B_v/B_z = 0.6 and decay more gradually in the lower B_v/B_z cases. SOL measurements typically show that the ratio T_i/T_e increases with the radius.⁴⁹ We see this trend in Fig. 8 for B_v/B_z = 0.3 and 0.2 but not for B_v/B_z = 0.6. This reversed trend for B_v/B_z = 0.6 is likely connected to the relatively flat T_e at large x. In the B_v/B_z = 0.6 case, the low-amplitude source of ∼33 eV electrons at large x [see the form of the plasma source, Eq. (5)] could be setting T_e in this region.

The flat T_e at large x could also be an artifact from the electrons running into a floor in the temperature at large x. However, we note that the minimum electron temperature allowed on our present grid is $Te,min=(2/3)T⊥e,min+(1/3)T∥e,min=11 eV$ based on $T⊥e,min=16 eV$ and $T∥e,min=1.1 eV$⁠, which is somewhat lower than T_e seen in this region. We can test this in future work by running higher-resolution runs, including a variable μ grid to better resolve low energies or by using exponential reconstructions, which is currently being added to the code.

The normalized root-mean-square (r.m.s.) electron-density fluctuation level in the x–z plane is shown in Fig. 9. For all three values of B_v/B_z, the density fluctuation levels are the largest in the source region $|z|<Lz/4$⁠. The normalized density fluctuation levels in the B_v/B_z = 0.6 case are fairly uniform along the field lines, while they tend to fall off by about a factor of 2–3 towards the sheaths in the smaller B_v/B_z cases. This effect could be a result of the stronger influence of the sheath on the potential as the distance from the source to the sheath is decreased. The instantaneous snapshots of $ñe$ (not shown) indicate a strong $k∥=0$ component for the largest B_v/B_z cases, while more parallel structure is apparent in the smaller B_v/B_z cases.

The fluctuation statistics can also give information about the strength of the electron adiabatic response for each simulation. By assuming that the electrons are isothermal along field lines, the parallel force balance satisfies

where ϕ_sh and n_sh are the electrostatic potential and electron density evaluated at the sheath entrances and ϕ_mid and n_mid are the same quantities but evaluated at the midplane (z = 0 m). To compute the cross-coherence diagnostic,^53,79,80 ordered pairs $(eϕmid/Te,eϕsh/Te+ln (nmid/nsh))$ falling in the region 1.318 m ≤ x ≤ 1.326 m (approximately where the maximum density and potential fluctuations are) are sampled at 1 μs intervals over a ∼1 ms period for each simulation. Figure 10 shows the resulting plots (normalized bivariate histograms), which all indicate a strong correlation between the two sides of Eq. (9), and so, the electrons are strongly adiabatic, meaning that the electron distribution function along a field line closely follows a Boltzmann distribution.⁸¹ This finding indicates that it might be possible to obtain similar results using a two-dimensional turbulence model (with reduced parallel dynamics and sheath-model boundary conditions) for the parameters considered here. To quantify the degree of non-adiabaticity, we define the parameter

where $ϕ¯mid$ is the left-hand side of Eq. (9) and $ϕ¯ad$ is the right-hand side of Eq. (9). The ϵ_na parameter measures the fraction of fluctuations in ϕ which are due to non-adiabatic effects. We find that ϵ_na is 0.094 for B_v/B_z = 0.6, 0.226 for B_v/B_z = 0.3, and 0.310 for B_v/B_z = 0.2, which is consistent with our expectation that the electrons become less adiabatic as B_v/B_z is decreased.

Figure 11(a) shows the radial profile of the autocorrelation time τ_ac (computed from time traces of the density fluctuations). In the SOL of the simulation, τ_ac tends to increase with the radius, which is a trend observed in measurements on NSTX (see Fig. 12 of Zweben et al.⁷⁰). The autocorrelation time for the B_v/B_z = 0.2 and B_v/B_z = 0.3 cases is found to vary between ∼5 μs and ∼9 μs, while the autocorrelation time for the B_v/B_z = 0.6 case exhibits a larger variation in the SOL, with τ_ac ≈ 4 μs for x < 1.34 m and increasing to ≈12 μs at the outer radial boundary. The autocorrelation times we observe in our simulations are lower than the τ_ac ∼ 10–40 μs reported by Zweben et al.⁷⁰ for the NSTX edge and SOL but are well within the τ_ac ∼2–20 μs range that is typical for edge and SOL turbulence in other tokamaks.^28,48

Figure 11(b) shows the poloidal and radial correlation lengths (L_pol and L_rad, respectively) using the electron-density fluctuations near the z = 0 m plane. The correlation length at a radial location is obtained by averaging the correlation length computed at several points in y. At an individual point, the correlation length is determined from the correlation function, which is constructed by computing the equal-time two-point autocorrelation function for density fluctuations separated by some distance Δy for L_pol or Δx for L_rad. Having observed a significant wave feature in the poloidal correlation function, we determined L_pol by fitting the poloidal correlation function to $e−|Δy|/Lpol cos (kwaveΔy)$⁠. The radial correlation function, which does not have a wave feature, is computed using the full width at half maximum (FWHM) as L_rad = FWHM/(2 ln 2).

For all three values of B_v/B_z, we observe that the ratio L_pol/L_rad is between 1.2 and 1.6 for most of the radial domain, which is similar to L_pol/L_rad ∼ 1–2 that is typically observed in tokamaks and stellarators.^28,48 An average L_pol/L_rad = 1.5 ± 0.1 was reported for representative Ohmic NSTX discharges,⁷¹ although larger ratios L_pol/L_rad ∼ 3–4 have been observed in some experiments⁸² and simulations.¹² 

There are two kinds of sheath-model boundary conditions that are commonly used in fluid and gyrokinetic codes. Logical-sheath boundary conditions enforce j_∥ = 0 at the sheath entrances, while current fluctuations into the sheath are permitted in conducting-sheath boundary conditions. Figure 12 shows the radial profiles of the steady-state parallel current into the sheath for the three cases under consideration. The currents have been normalized to peak steady-state ion saturation current j_sat = q_in_ic_s, where $cs=(Te+γTi)/mi$ and γ = 3 is used because the collisionless layer in front of the sheaths should be resolved in all three cases. All three cases are quite quantitatively similar, and the outward sheath currents are found to be highly symmetric in z, which is consistent with the strong adiabatic response shown in Fig. 10. A large excess electron outflow (negative current) is seen in the hot source region (near x = 1.3 m), which is compensated by a large excess ion outflow (positive current) just outside the source region. The peak values are approximately 20% of the ion saturation current, which motivates future studies regarding how the use of various sheath-model boundary conditions affect turbulence in these simulations.

We have developed a model to investigate curvature-driven SOL turbulence in a simplified helical-magnetic-field geometry and performed numerical simulations of the system using an electrostatic gyrokinetic continuum code. The blobs in our simulations appear to originate as radially elongated structures that extend from the source region into SOL and get broken up by sheared poloidal flows. The blobs appear to efficiently transport plasma across the magnetic field, leading to radial particle fluxes that are much higher than Bohm-flux estimates. Such large-amplitude and large-scale blobs were not observed in a set of simulations we performed without magnetic-curvature effects. We note, however, that coherent structures with high plasma density have been observed in linear devices with the negligible magnetic curvature.^83,84 The mechanism that polarizes such coherent structures in linear devices and leads to outward radial propagation could be due to neutral wind.⁸⁵ 

We characterized the turbulence using a variety of diagnostics and found that various quantities of interest are within the range expected for SOL turbulence in tokamaks, such as fluctuation levels, autocorrelation times, and correlation lengths. A summary of some quantities observed in our simulations is given in Table I, which also includes experimental values from NSTX SOL.^24,70 We know that there are a number of important physical effects (e.g., complete magnetic geometry, magnetic fluctuations, and atomic physics) that need to be added to the simulations in order to expect quantitative accuracy for detailed comparisons with experiments, but it is interesting to see that the present simulations are already in the right ballpark qualitatively. This and other recent work indicate the general feasibility of using continuum codes to simulate gyrokinetic turbulence in the edge and SOL regions of tokamaks.

Even in this simple limit, we began to explore a number of physical processes. We varied the magnetic-field-line pitch in a set of simulations, which indicated an increasing level of radial turbulent particle transport with the decreasing pitch. A cross-coherence diagnostic comparing potential fluctuations at the sheaths with those at the midplane indicated that all three simulations appeared to fall into a similar turbulent regime with strongly adiabatic electrons. The application of this model to investigate turbulence in the Helimak device^44,47 has also been performed and will be reported elsewhere.

The helical-SOL model can be extended by the addition of a closed-magnetic-field-line region (with periodic boundary conditions in the parallel direction). While the Gkeyll code can already perform simulations with periodicity in the parallel direction, additional work is required to simultaneously include both open and closed-magnetic-field-line regions in the same simulation. The addition of good-magnetic-curvature regions and electromagnetic effects are also important extensions that will make this model more applicable to tokamaks. Since our model is relatively simple compared to a realistic tokamak SOL, the helical-SOL model could also eventually serve as a test case for the cross verification of gyrokinetic boundary-plasma codes. This test case might be useful for revealing major discrepancies due to different numerical approaches, sheath-model boundary conditions, and collision operators implemented in various codes relatively early on in the development cycle before more significant investments are made.

See supplementary material for a video comparing the late-time electron-density (in 10¹⁸ m⁻³) dynamics in the plane perpendicular to the field lines at z = 0 m for the B_v/B_z = 0.6 and B_v/B_z = 0.3 cases.

We thank S. Zweben for useful discussions about NSTX SOL measurements and J. Juno for setting Gkeyll up on the Stampede cluster. E.L.S. would also like to acknowledge useful discussions with J. Nichols concerning ion temperature measurements in the SOL and S. Zweben and M. Kunz for providing feedback on the manuscript. This work was funded by the U.S. Department of Energy under Contract DE-AC02-09CH11466, through the Max-Planck/Princeton Center for Plasma Physics and the Princeton Plasma Physics Laboratory. E.L.S. prepared this manuscript in part under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. G.W.H. and A.H. were supported in part by the SciDAC Partnership for Multiscale Gyrokinetic Turbulence. A.H. was also supported in part by the Laboratory Directed Research and Development program. Some simulations reported in this paper 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. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation under Grant No. ACI-1548562.

This document was prepared as an account of work sponsored by an agency of the U.S. government. Neither the U.S. government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.

We consider a problem in which a uniform mass source S_ρ and energy source S_E are continuously active in the region $|z|<LS/2$⁠. This fluid flows out to perfectly absorbing boundaries at $|z|=Lz/2$⁠. We treat the plasma as a single fluid with mass density ρ ≈ n_em_i, pressure p = n_e(T_e + T_i) = 2n_eT (where T is an average of the electron and ion temperatures), and energy density (3/2)n_e(T_e + T_i), so S_ρ = m_iS_n and S_E = 3T_srcS_n, where S_n is the electron and ion particle source rate and T_src = (T_e,src + T_i,src)/2 is the effective single-fluid source temperature. This system is described by the steady-state ideal fluid equations (neglecting thermal conduction and viscosity)

where u is the fluid velocity, ρ is the mass density, and p is the pressure. We treat the source as having no mean flow in the z direction.

We integrate these equations from z = 0 to an arbitrary position z < L_S/2 and use the boundary condition u(z = 0) = 0 to get

where p₀ ≡ p(z = 0). The first two equations can be solved for ρ and p respectively, and we obtain a quadratic equation for u(z) by substituting these expressions into the last equation. The solution to this system is

Since the pressure cannot be negative, the only physical solution for small z is the negative branch for u(z) and the positive branch for p(z). The central pressure p₀ is determined by the boundary conditions at $|z|=Lz/2$⁠. A steady-state solution at a perfectly absorbing wall requires $M≥1$ at the wall,⁸⁶ where the Mach number $M(z)≡u(z)/cs(z)=u(z)/(5/3)p(z)/ρ(z)$⁠. The $M≥1$ requirement is equivalent to the Bohm criterion for a steady-state sheath. We see that

The maximum possible value of $M(z)$ occurs at z that minimizes p(z). This value z_max turns out to be z that makes the radicand in Eq. (A7) zero, so we find that

This says that the largest possible value of $M$ is 1 [when z has been made as large as possible, as given by Eq. (A11)]. This (barely) satisfies the outflow requirement that $M≥1$ at a perfectly absorbing wall. Equation (A11) then provides a constraint on the value of p₀ such that $M=1$ is achieved at the end of the source region, z_max = L_S/2

Using this expression for p₀, we have the following profiles in the source region $0<|z|<LS/2$

In order to use these profiles to initialize a Maxwellian initial condition for a kinetic simulation, we note that these profiles correspond to density (n = ρ/m_i) and temperature (T = m_ip/(2ρ)) profiles in the source region given by

where $css=(5/3)2Tsrc/mi$⁠. In the source-free regions z > L_S/2 or z < L_S/2, n(z), T(z), and u(z) are all constant and equal to the value that their respective profiles evaluated at the corresponding edge of the source region at z = L_S/2 or z = –L_S/2. The 1-D equilibrium profiles Eqs. (A16), (A18), and (A19), the density source in the helical-SOL simulations Eq. (5), and the temperature profiles of the electron and ion sources are used to generate spatially varying initial conditions in (x, y, z).

One could go further by calculating the slight difference between the ion-guiding-center-density and electron-density profiles which gives the desired equilibrium potential ϕ(x, y, z) when the gyrokinetic Poisson equation is solved. For now, we simply set $nig(x,y,z)=ne(x,y,z)$ and initialize with ϕ = 0, as we did in the LAPD simulations of Shi et al.¹⁴