Gyrokinetic simulation of driftwave instability in field-reversed configuration

A field-reversed configuration (FRC) is a compact toroid with magnetic field predominantly in the poloidal direction. Somewhat after the discovery of FRCs in the late 1950s,^1,2 it was recognized that their characteristic high-β (ratio of plasma pressure to magnetic pressure) and compact size would be favorable to the design of a fusion reactor. Despite significant advances, maintaining macrostability remained a challenge throughout resurgence in research that took place in the 1970s and 1980s.³ Renewed interest in FRC research in the late 1990s onward leads to the development of new formation techniques and better theoretical understanding of macroscopic instabilities.⁴ 

The “advanced beam-driven” flavor of the FRC was originally conceived in 1993 by Rostoker et al.^5,6 The salient feature of this concept is a dominant nonadiabatic fast ion population, which brings a number of stability benefits to the configuration. It was known at the time that nonadiabatic ions in tokamaks were relatively insensitive to microturbulence,⁷ likely due to the spatial averaging effect of ion Larmour radii comparable to the plasma size. It was hypothesized that a significant fast ion population, in an FRC, could also stiffen the plasma against some macroinstabilities (e.g., tilt mode). Other suppression techniques, such as electric biasing schemes or additional magnetic windings, could be used to stabilize any remaining macroscopic modes. In practice, neutral beam injection (NBI) is used to provide the fast ion population, which is the origin of the “advanced beam-driven” moniker.

The fast-ion conjecture was taken up by Tri Alpha Energy, Inc. (TAE), which began work on FRC experiments in the late 1990s. The C-2 experiment at TAE began in 2008, and among other goals, aimed to demonstrate that NBI could stabilize and sustain an FRC plasma.⁸ From 2008 to the present, remarkable progress in FRC plasma confinement in C-2 and its successor, C-2U, has brought the advanced beam-driven FRC concept to a point where microturbulent instability is one of the foremost physics obstacles.^9–11 

In this paper, we begin to investigate microturbulence in advanced beam-driven FRC via gyrokinetic simulations. Section II summarizes the state of recent progress at TAE and motivates kinetic turbulence simulation. The initial simulation model is electrostatic and linear. Details of this simulation model are presented in Section III. Remarkably, ion-scale driftwave turbulence is found to be linearly stable in the FRC core. These results and others are presented in Section IV. Interpretation of results, possible stabilization mechanisms, and caveats are discussed in Section V. This work represents a first effort to characterize microturbulence in realistic FRC geometry, but model refinements and additional physics are needed. Future work and direction of code development are discussed in Section VI.

The C-2 experiment at TAE was designed as a testbed for advanced beam-driven FRC physics, with an ambitious list of goals including exploring effects of NBI, rapid repetition of discharges (≤10 min), reproducible plasmas with lifetimes of at least 5 ms (longer than any characteristic decay times), and reproducibly macro-stable plasmas to allow the study of turbulent fluctuations and transport.⁸ From the first data taken on C-2 in 2008, to the end of the campaign in early 2014, dramatic strides were taken in terms of operational expertise, plasma quality, and lifetime. The culmination was the High Performance FRC regime (HPF) which combined synergies of NBI, electric end-biasing on open field lines, and advanced wall cleaning techniques to achieve reproducible plasma discharges with lifetimes exceeding 5 ms.¹¹ 

The upgraded experiment, C-2U, which began its campaign in March of 2015, made a number of diagnostic and control upgrades to C-2. The most critical of these extended the success of the HPF regime by upgrading the neutral beam system to six new 15 keV beams with a total power of 10 MW, significantly exceeding the 4 MW total beam power available on the original C-2 machine. In June of 2015, C-2U achieved sustainment of a plasma for over 5 ms, limited only by the neutral beam pulse length.¹² 

Success on the C-2 and C-2U campaigns in achieving macro-stable, sustained plasma has allowed first measurements of turbulent fluctuations to be taken.¹³ In the current physics regime (1 keV ion temperature, 10¹⁹cm⁻³ density), particle and energy confinement times scale favorably as temperature is increased.¹¹ The existence of an upper temperature limit to this favorable scaling remains an open question, which is now possible to address due to recent substantive advances in macro-stability control.

Experimental studies of turbulent transport in FRCs have previously been limited by short experimental confinement times. Nonetheless, significant investigations of FRC transport characteristics have been carried out. In most experimental schemes, the FRC is isolated from material walls, and the particle recycling rate is low. To date, beam fueling in real experiments has been absent or small compared to the bulk ion population. Radiative losses and work from magnetic compression are also typically small. Consequently, the particle and energy confinement times are well approximated by the global particle and energy decay times, respectively. It has also been observed, despite modest collisionality, that electron thermal transport is dominated by convection rather than conduction.¹⁴ Because energy is primarily conveyed by particle convection, the energy confinement is governed by the particle confinement.

The nature of confinement varies qualitatively in different regions of the plasma. On the closed field lines of the FRC core, confinement is determined by cross-field transport. In the 1970s, the lower hybrid drift instability (LHDI) was identified as the most linearly unstable electrostatic micro-instability in this region; however, early numerical simulations of LHDI^15–17 overpredict experimentally observed density fluctuations. Electromagnetic micro-instabilities may play a larger role in cross-field transport. Because the equilibrium magnetic field approaches zero at the O-point, even small magnetic perturbations can connect flux surfaces, allowing electron convection. The need to understand the transport contributions of electrons, in detail, motivates inclusion of electron kinetic responses in the physics simulation model.

Outside the core, it is critical to understand transport across the separatrix, which acts as a boundary condition for both the core and open field line scrape-off layer (SOL). Significant plasma rotation can create a complex separatrix structure not captured by single-fluid MHD, where each species, α, has unique characteristic surfaces, P_α,¹⁸ 

Minimally, a two-fluid or kinetic model is required to capture this behavior.

Confinement in the SOL is qualitatively very different from the core and separatrix. Particle transport is a balance of radial diffusion and parallel flow on the open field lines. Parallel flow is strongly influenced by end conditions, which may include field constrictions (e.g., magnetic mirrors), plasma sheaths, and end biasing schemes at the divertor wall. Plasma in the SOL region is not cold. Observations indicate that the temperature in the SOL, just outside the separatrix, is comparable to the core, indicating strong coupling between the regions. As a result, the simulation domain in any transport calculation should include core, separatrix, and SOL.

The challenge of quantitatively predicting confinement scaling in an advanced beam-driven FRC will require modern high performance computing (HPC). To resolve driftwave scale turbulence in a domain including the coupled FRC core and SOL regions, a predictive simulation must capture length scales from the device size, ∼1 m diameter and ∼10 m in length, down to the length scale of the density gradient, which may be as small as 6 × 10⁻⁵ m just outside the separatrix. Even finer grid resolution may be required. Recent measurements of density fluctuations indicate that electron-scale turbulence, with wavelengths significantly shorter than the ion Larmour radius, dominates in the FRC core.¹³ Relevant timescales for the problem range from mesoscale MHD modes down to the electron transit time on a field line in the FRC core. The length of the C-2U device divided by the sound speed is on the order of ∼100 μs, while the electron transit time in the core is ∼0.2 μs. The need to capture many physics features, including large ion orbits, complex separatrix structure, and electron contributions to both energy transport and turbulence mitigation, motivates use of a kinetic model. Characterizing such a wide range of spatial and temporal scales represents a significant computational challenge, which, to date, is only being addressed by HPC. In the following work, we begin efforts to address this challenge by applying a highly scalable, mature particle-in-cell (PIC) code, the Gyrokinetic Toroidal Code (GTC), to the problem of FRC microinstability.

To the authors' knowledge, first-principles simulation of turbulent transport using realistic FRC magnetic geometry has never been previously carried out. In this study, we apply the Gyrokinetic Toroidal Code (GTC)¹⁹ to driftwave turbulence using numerical geometry based on the TAE advanced beam-driven FRC experiment, C-2. GTC is a full-featured, well-benchmarked, kinetic PIC code with gyrokinetic or fully kinetic ion models, kinetic electron effects,²⁰ electrostatic or fully electromagnetic field solvers, a Fokker–Planck collision operator,²¹ perturbative δf or full-f resolution of the particle distribution functions, and an efficient, fully 3D magnetic coordinate grid which follows either analytically or numerically specified magnetic equilibrium geometry. GTC has been applied to simulations of microturbulence,²² energetic particle transport,²³ Alfvén eigenmodes,^24,25 and to MHD modes include tearing²⁶ and kink²⁷ instabilities in tokamak plasmas. A number of upgrades were recently implemented in GTC for the purpose of FRC simulation. These include (1) a coordinate transformation pipeline to produce Boozer magnetic coordinates from numerical FRC geometry, and (2) reimplementation of the electrostatic field-solver in order to optimize parallelism for an equilibrium magnetic field line orientation in the poloidal direction.²⁸ 

For the initial effort presented in this paper, only linear mode growth of electrostatic instabilities is considered. Electromagnetic turbulence may play a crucial role in transport and will be investigated in the future. The ion model is gyrokinetic and electrons are drift kinetic using a realistic mass ratio. Boundary conditions for the particles are reflective on the radial boundaries and periodic on the axial boundary. The field solver assumes a fixed value field on the edges of the grid. Simulations are initialized with particles distributed uniformly across computational cells, with random spatial placement within each cell. The velocity distribution of particles is assumed to be Gaussian about the thermal temperature for that species. Both fields and particle distributions are computed with the perturbative δf model. Early simulations indicate the effect of collisions on linear growth rate and mode frequency is small, consistent with a low effective collisionality (see Section IV). For simplicity, collisions are excluded in the following results, but are under investigation and will be reincorporated in future simulations.

The magnetic coordinate system in GTC provides numerical advantages in the required grid size for the field solver, however is unable to simulate the separatrix, where there is a coordinate singularity. To circumvent this problem, the core and SOL regions have been simulated separately with the separatrix excluded. The magnetic axis is also excluded, since the gyrokinetic model is invalid at the field null in the center. A further approximation is that, in linear simulations, the toroidal wavelength is much shorter than the radial wavelength. It may be expressed in terms of the toroidal and radial wavenumbers as k_r ≪ k_ζ, and allows single flux surfaces to be simulated in isolation. The result is significant savings in simulation time due to reduced domain and simplified Laplacian operator. The approximation is accomplished by modifying particle gather and scatter operations to only sample one flux surface radially. Both core and SOL regions, with the simulated flux surfaces indicated, are shown in Fig. 1.

While the end goal of this work is to produce global domain simulations of an FRC, using a significantly smaller domain for initial simulations results in a rapid turnaround time for individual simulations. This rapid turnaround is critical to enable careful evaluation of all aspects of the physics model before more-expensive global simulations are carried out. Other approximations made in this paper are motivated similarly. Rapid iteration of simulations is an absolute necessity in the development of a complete integrated physics model for FRC transport.

In this section, we report on initial results of driftwave instability simulations in FRC geometry. The simulation in the FRC core is performed on a flux surface with $r/a≈0.906$⁠. In the SOL, a flux surface with $r/a≈2.244$ is used. The radial coordinate, r, is the position of the flux surface, measured outward on the midplane from the magnetic axis, and a, the minor radius of the FRC, is the value of r at the magnetic separatrix. These locations were chosen by performing radial scans to determine flux surfaces with maximal growth rates.

Driftwave instabilities are driven by density and temperature gradients. These gradients may be characterized by their length scales, where, for a plasma parameter, g, the inverse length scale is defined as $1/Lg≡∂∂rln(g)$⁠. Both temperature and density gradients are included. The ratio of the scale lengths of these two gradients, η, indicates their relative importance. For ions and electrons, $ηi=LTi/Ln$ and $ηe=LTe/Ln$⁠, respectively. In all of the following simulations $η≡ηi=ηe=1$⁠. The simulated ion species is deuterium. The major radius, measured from the machine axis to the magnetic axis, is $R0=26.96 cm$⁠, and the minor radius is $a=11.15 cm$⁠. Simulation parameters and derived quantities are presented in Table I, including the ion and electron gyro-radii (⁠ $ρi=miTi/(eB), ρe=meTe/(eB)$⁠).

Transit frequencies and effective collisionalities are listed in Table II. The transit frequency describes how often a particle completes a circuit on a field line. For electrons and ions on a field-line of length l, the respective transit frequencies are defined $νtr,e=Vth,e/l$ and $νtr,i=Vth,i/l$⁠, respectively, where $Vth,α$ is the thermal velocity of a particle species, α. Transit frequencies are normalized by the ratio of the major radius, R₀ to the sound speed, $Cs=Te/mi$⁠. In the core, the field-line length is estimated by $l=(R0+r)×π$⁠, and in the SOL, $l≈5 m$⁠, which is the length of the simulation domain. Effective collisionality is the ratio of the collisional frequency to the transit frequency and characterizes the relative importance of collisional effects. When $να1−α2*≲1$⁠, collisions have only a small effect and become very important as $να1−α2*≫1$⁠. Here $να1−α2*$ is the effective collisionality between species α₁ and α₂, and the unstarred $να1−α2$ represents the unnormalized collisional frequency. For electrons, e, and a single ion species, i, $νe−e*=νe−e/νtr,e, νe−i*=νe−i/νtr,e$⁠, and $νi−i*=νi−i/νtr,i$⁠. Notably, the effective collisionality in an ideal FRC is quite low compared to a tokamak with the same major radius and temperature. Because the FRC has equilibrium magnetic field only in the poloidal direction, the field line length is short compared to a tokamak where field lines must make many toroidal circuits before closing. As a direct result, the transit frequency in the FRC is much higher, and, in turn, the effective collisionality is much lower.

In the core, ion-scale driftwave modes are found to be strongly stable using experimentally realistic values for density and temperature gradients. Artificial increase of these gradients to the limits of gyrokinetic validity still yields no instability in the core. Magnetic field variation on the core-simulation flux surface is shown in Fig. 2. Modes with higher toroidal wavenumber, $kζρi$⁠, well into electron-scale lengths, are easier to excite in the simulation, but these modes do not typically contribute significantly to transport. Careful investigation of this high-wavenumber electron-scale driftwave turbulence will be reported in a future paper.

In the scrape-off layer, realistic values of the ion temperature gradient inverse length scale are $R0/LTi≈2.5$⁠. As shown in Fig. 3, modes with finite growth rate are found in the SOL using temperature and density gradient drives in the range $R0/L≈2−6$⁠. By varying the instability drive via $R0/L$⁠, it is possible to detect the minimal drive for which a linear mode becomes unstable. This linear instability threshold is comparable (with some caveats) to critical gradient instability thresholds observed in the experiment. For ion-scale length modes, the linear stability threshold occurs for a drive strength in the range of $R0/L≈3−4$⁠, which is qualitatively consistent with nonlinear critical gradient measurements, indicating $R0/Lcrit≈5−6$⁠, taken in recent Doppler backscattering data on the C-2 experiment.¹³ Turbulence into electron-scale lengths, represented by green curves to the left in Fig. 3, is also present in the SOL. As expected, these higher wavenumber electron-scale length modes are easier to excite and, correspondingly, have lower linear stability thresholds.

The dispersion relation for a typical SOL mode, with a driving gradient $R0/L≈4$⁠, is shown in Fig. 4. Notably, the sign of the mode frequency is consistent with both the electron diamagnetic frequency as well as the ion grad-B frequency. As expected, the growth rate trends upwards as the gradient strength, $R0/L$⁠, is increased.

The corresponding mode structure for the SOL instability, for an intermediate scale length mode (⁠ $kζρs=5.48$⁠), is shown in Fig. 5. The ρ_s length normalization is similar to ρ_i but includes electron temperature as well as ion temperature and is defined $ρs=mi(Ti+Te)/(eB)$⁠. The mode amplitude in the top panel displays mixed $k∥=0$ and low finite- $k∥$ contributions. Although the field line has bad-curvature everywhere (second panel), the length scale of the magnetic field gradient (fourth panel) is comparable or larger everywhere in the simulation domain, so the drift terms have an overall stabilizing effect.

In this paper, electrostatic gyrokinetic simulations of driftwave turbulence in FRC geometry have been performed. Notably, for ion scale lengths (⁠ $kζρi≲1)$⁠, no instabilities exist in simulation in the FRC core region. A number of mechanisms may contribute to this observed stability. These candidates for core stabilization are under investigation and include finite Larmour radius, magnetic field gradient, and other effects of the magnetic geometry.

Driftwave turbulence is still marginal in the SOL, but possible to drive with realistic temperature and density gradients. A typical dispersion relation in the SOL indicates increased growth with stronger drive and indicates electron diamagnetic response or ion grad-B response as energy channels for instability. Further analysis of particle phase space resonances is underway and will help complete the physical picture of the instabilities in the SOL.²⁹ Realistically, modelled input equilibrium suggests that the peak driving gradients occur in the SOL region, just outside the magnetic separatrix.

These simulation results are consistent with the recent experimental data taken from the C-2 experiment.¹³ Doppler backscattering measurements of density fluctuations indicate a unique inverted turbulence spectrum with ion-scale lengths inside the separatrix. Measured fluctuation amplitude outside the separatrix is significantly higher, although still quiescent in the C-2 HPF operating regime, and displays a more-typical exponential turbulence spectrum.

Despite this consistency, GTC simulation results are not directly comparable to turbulent fluctuation measurements. Linear growth rates in simulation do not translate directly to nonlinear fluctuation amplitudes. Different nonlinear saturation mechanisms may come into play for separate linear modes, resulting in saturation levels independent of growth rate. A fully nonlinear simulation is necessary to capture these saturation mechanisms and directly compare simulation to experiment. These and other limitations of the simulation model are caveats to the results presented here, and provide a roadmap for future investigation and code development.

Towards the goal of predictive transport simulations of an FRC-based fusion reactor, a number of additional physics features need to be included in a working simulation model. One of the major challenges is to address the nature of coupling between the core, separatrix, and SOL regions of the advanced beam-driven FRC. The magnetic coordinates used in GTC and in other highly optimized PIC codes are singular at the magnetic separatrix, making simulation at this location numerically intractable. To address the separatrix problem, development of A New Code (ANC), which leverages GTC physics but uses a cylindrical coordinate grid, is well underway. Currently, the ANC model is perturbative (δf) and includes a drift kinetic ion species and adiabatic electrons, which are parallelized in groups. A fully 3D field solver is parallelized and spectrally decomposed in the azimuthal coordinates. The same realistic numerical equilibria used to produce magnetic equilibria in GTC are read directly into ANC. The cylindrical coordinate grid sacrifices some efficiency, but spans the critical separatrix region.

Other mandatory physics features (see Section II) include electromagnetic turbulence and fully kinetic ions. These simulation capabilities already exist in GTC using the magnetic coordinate grid. Development efforts will focus on porting GTC's proven physics algorithms into ANC. Simultaneously, physics analysis can move forward using existing GTC features but restricted to simulation domains that do not cross the separatrix.

Beyond inclusion of fully kinetic ions, a cross-separatrix simulation domain, and electromagnetic perturbations, long term development goals include realistic end boundary conditions on open field lines, and the use of more general FRC equilibria which allow finite toroidal magnetic field. Code features from the new code, ANC, are also intended to be integrated into the GTC suite to produce a more robust, integrated physics code.

As mentioned in Section V, fully nonlinear simulations are required to quantitatively predict transport as well as to make direct comparison between simulation and experiment. The physics model used in both GTC and the newly developed ANC is fundamentally nonlinear. Nonlinear physics will be investigated after the spectrum of linear modes is more thoroughly understood.

This work was carried out as part of the author's thesis work at University of California, Irvine, and continued at Tri Alpha Energy, Inc. Work at University of California was performed with the support of the Norman Rostoker Fellowship (Grant No. TAE-200441) and by the DOE SciDAC GSEP center. GTC simulations used resources on the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory (DOE Contract No. DE-AC05-00OR22725), and the National Energy Research Scientific Computing Center (DOE Contract No. DE-AC02-05CH11231).