11in

# Zonal Flow Excitation in Electron-Scale Tokamak Turbulence

Stefan Tirkas  0000-0002-8050-617X<EMAIL_ADDRESS>University of
Colorado at Boulder, Boulder, CO 80309, United States of America  Haotian Chen
0000-0003-3896-4369 Southwestern Institute of Physics, Chengdu 610041, China
Gabriele Merlo  0000-0003-4877-1456 The University of Texas at Austin, Austin,
TX 78712, United States of America  Frank Jenko  Max Planck Institute for
Plasma Physics, Boltzmannstr 2, 85748 Garching, Germany  Scott Parker
Renewable and Sustainable Energy Institute, University of Colorado at Boulder,
Boulder, CO 80303, United States of America

###### Abstract

The derivation of an intermediate-scale gyrokinetic-electron theory in
nonuniform tokamak plasmas [Chen H. et al 2021 Nucl. Fusion 61 066017] has
shown that a Navier-Stokes type nonlinearity couples electron-temperature-
gradient (ETG) modes and zonal flow (ZF) modes with wavelengths much shorter
than the ion gyroradius but much longer than the electron gyroradius. This
intermediate-scale ETG-ZF coupling is typically stronger than the Hasegawa-
Mima type nonlinearity characteristic of the fluid approximation and is
predicted to lead to relevant zonal flow generation and ETG mode regulation.
Electron-scale, continuum, gyrokinetic simulation results are presented here
which include both single-mode ETG and full-spectrum ETG turbulence. The zonal
flow generation due to single ETG modes is investigated and the single-mode
intermediate-scale results are found to be in agreement with theory. The full-
spectrum results are then presented and explained qualitatively in terms of
the single-mode results. It is found that the ETG-driven zonal flows regulate
intermediate-scale electron heat flux transport to levels in the predicted
range.

##  I Introduction

Electron-scale turbulence is a plausible explanation for the anomalous
transport of electron energy well above the neoclassical values seen in a
variety of tokamak plasma scenarios BatchTrans ; ConfineITER ; HortonTransport
; DoyleTransport . Additionally, electron energy transport may become more
important in future burning plasma experiments such as ITER because the
electron channel is preferentially heated by Coulomb collisions with fusion
alpha particles. The electron-temperature-gradient (ETG) instability produces
radially-elongated streamers at the electron gyroradius scale and is a primary
candidate to explain electron-scale transport ConfineITER ; DoyleTransport ;
LeeETG ; DorlandJenkoETG1 ; DorlandJenkoETG2 ; DorlandJenkoETG3 ; Nevins06 .
Electron heat flux due to ETG turbulence has been seen to play a role in
various tokamak experiments RenETG ; GriersonETG ; RyterETG ; KieferETG , and
the inclusion of electron-scale dynamics at the ion scale has resulted in
better agreement with experimental heat flux levels HowardLmode ; HollandHmode
.

Electron-scale ParkerZF ; ColyerETG  and multiscale HowardLmode ; HollandHmode
long-time, large-box gyrokinetic flux-tube simulations have reported that
intermediate-scale zonal flows (ZFs) help to regulate streamer turbulence in
the quasi-saturated state and can eventually become dominant. These results
are inconsistent with fluid ETG turbulence models in which zonal flow
generation occurs via the standard Hasegawa-Mima nonlinear mechanism
HasegawaMima , which is significantly weaker for ETG turbulence than for ion-
temperature-gradient (ITG) turbulence LinFluidETG ; ChenFluidETG . Moreover,
while shearing due to zonal flows generated by ITG turbulence can suppress ion
transport levels, the finer scale ETG turbulence is unlikely to be affected by
the ITG-driven zonal flows BatchTrans ; KimETG . These effects have led to the
expectation of a streamer-dominated steady state at the electron scale.

A weak-turbulence, toroidal, gyrokinetic-electron analysis Chen_ETG_ZF  of
nonuniform tokamak plasmas in the intermediate-scale
($k_{\perp}^{2}\rho_{e}^{2}\ll 1\ll k_{\perp}^{2}\rho_{i}^{2}$) results in a
Navier-Stokes type nonlinearity which is typically stronger than the Hasegawa-
Mima coupling of the fluid approximation. Here $k_{\perp}$ is the wavenumber
perpendicular to the magnetic field and $\rho$ is the particle gyroradius. The
subscripts $i$ and $e$ denote ions and electrons respectively. This stronger
Navier-Stokes type nonlinearity is predicted to drive notable ZF generation
and ETG mode regulation at intermediate scales when compared to the Hasegawa-
Mima type coupling of the short-wavelength fluid regime Chen_ETG_ZF .

ETG mode saturation in the short-wavelength fluid regime has been extensively
studied and includes saturation mechanisms such as secondary instabilities
DorlandJenkoETG1 ; DorlandJenkoETG2  and toroidal inverse-cascading
LinFluidETG ; ChenFluidETG ; KimETG  which would lead to a turbulent state
characterized by intermediate-scale ETG modes. After the initial transition to
the intermediate scale, the nonlinear interaction between ETG and ZF modes is
expected to be enhanced such that zonal flows may grow to regulate long-term
steady state transport levels as measured by experiment Chen_ETG_ZF . The
intermediate-scale turbulent state might then be characterized by kinetic
saturation mechanisms such as standard quasilinear estimates JenkoSat1 ;
JenkoSat2 ; BourdSat  and $\textbf{E}\times\textbf{B}$ particle trapping
ParkerITG .

Figure 1: GENE linear simulation results showing growth rates and real mode
frequencies for (a) CBC ITG benchmark case, (b) electron-scale ETG case, and
(c) ion-scale ETG case showing CTEM modes as well. Positive(negative)
frequencies indicate propagation in the ion(electron)-diamagnetic direction.

Collisionless Cyclone base-case (CBC) simulations of ETG turbulence initially
reported the algebraic growth of zonal flows into late times ParkerZF .
Electron-scale MAST simulations ColyerETG  demonstrated that the long-time
saturated electron heat flux scales roughly proportionally to the
collisionality, and this was tied to the nonlinear interaction of ETG modes
with zonal flows which are well-known to be damped by collisions DiamondZonal
. In both cases, an initial turbulent state developed characterized by ETG
mode streamers which were eventually suppressed by the slow growth of
intermediate-scale zonal flows. DIII-D and Alcator C-Mod simulations
HollandHmode ; HowardLmode  involving multiscale ion and electron dynamics
also saw significant intermediate-scale zonal flow generation which helped to
suppress ITG and ETG turbulence into the late stage.

This paper provides results which compare the generation of zonal flow by ETG
turbulence in electron-scale, gyrokinetic simulations with theoretical
predictions in the intermediate-scale limit Chen_ETG_ZF . We first provide the
details of the simulation parameters in Section II. This is followed by
analysis of two types of nonlinear simulations. The “single-mode” results
serve to illuminate the role of a single ETG mode in generating zonal flow and
are covered in Section III. As the theoretical description is limited to a
single ETG mode for tractability Chen_ETG_ZF , these results convey the
primary scope of this paper. The “full-spectrum” nonlinear simulations
provided in Section IV include a typical range of ETG modes and are
qualitatively explained in terms of the single-mode results. In both types of
simulations it is found that intermediate-scale zonal flows are primarily
driven by slowly-saturating intermediate-scale ETG modes. The results are in
good agreement with the gyrokinetic theory in the intermediate-scale and the
electron fluid models at the short wavelength scale.

##  II Simulation Model and Parameters

We employ GENE DorlandJenkoETG1 ; GeneCode , an Eulerian 5-d gyrokinetic
continuum code, in the flux-tube limit appropriate for electron-scale
turbulence. Gyrokinetic ions and electrons are taken with standard Cyclone
base-case (CBC) parameters which are typical of H-mode core plasmas, but here
a simplified circular geometry is used CBC_1 ; CBC_2 ; CBC_3 ; GorlerBenchmark
. First, the linear ion-scale benchmark in Ref. GorlerBenchmark  was verified,
then the simulation was converted to the electron scale by reducing the
perpendicular box dimensions by a factor of $\sqrt{m_{i}/m_{e}}\sim 42$. Here,
$m_{e}$ is the electron mass, and $m_{i}$ is the ion mass which is taken to be
the proton mass, $m_{p}$. For the electron-scale case the ion temperature
gradient was set to zero to suppress long-wavelength ion turbulence and focus
on electron-scale physics.

The mode frequencies, $\omega$, and growth rates, $\gamma$, resulting from the
linear GENE simulations are shown as functions of $k_{y}\rho_{i}$ in Fig. 1
above. The frequencies are normalized to units of $R/c_{s}$, as listed with
other GENE normalizations in Table I. Here, $k_{y}$ is the binormal wavenumber
of the GENE coordinate system, $R$ is the tokamak major radius, and
$c_{s}=\sqrt{T_{i}/m_{e}}$ is the ion sound speed with $T_{i}$ the ion
temperature. The ITG benchmark case is shown in Fig. 1 alongside the electron-
scale ETG case in Fig. 1. Ion-scale ETG results are shown in Fig. 1, where the
collisionless trapped-electron mode (CTEM) CBC_TEM ; AdamCTEM ; ChenCTEM  is
included. One can see the ITG mode in the lower $k_{y}\rho_{i}$ range of Fig.
1 characterized by propagation in the $\omega>0$ ion diamagnetic drift
direction. The ETG mode becomes unstable at higher $k_{y}\rho_{i}$ where
$\omega$ crosses to the negative electron diamagnetic drift direction
HortonTransport . While the CTEM is expected to contribute to electron
transport BatchTrans ; JenkoTEM ; XiaoTEM ; ChenCTEM2 , Fig. 1 shows that it
is stable in the intermediate-scale range.

Figure 2: Original (small-box, no collisions) nonlinear simulation results
showing (a) time-marked electrostatic electron heat flux and electrostatic
potential contours for the (b) early nonlinear phase (green marker), and (c)
late zonal phase (red marker).

Figure 3: Original (small-box, no collisions) nonlinear simulation results
showing Fourier-space electrostatic potential contours for the (a) initial
time (blue marker), (b) early nonlinear phase (green marker), and (c) late
zonal phase (red marker). The markers correspond to the times marked in Fig.
2.

The reference values and radial profiles are taken as specified in Ref.
GorlerBenchmark  with reduced values of $m_{i}$ and $\beta=8\pi
n_{e}T_{e}/B^{2}$. Here $n_{e}$ and $T_{e}$ are the electron density and
temperature, and $B$ is the on-axis magnetic field. $\beta=10^{-5}$ was chosen
to keep the simulation nearly electrostatic in order to avoid transport due to
electromagnetic fluctuations. The safety factor, $q=rB_{\phi}/RB_{\theta}$,
and magnetic shear, $\hat{s}$, profiles are given by
$q(r)=2.52(r/a)^{2}-0.16(r/a)+0.86$ and $\hat{s}=\frac{r}{q}\frac{dq}{dr}$
GorlerBenchmark . $B_{\phi}$ and $B_{\theta}$ represent the toroidal and
poloidal magnetic field components respectively. The radial flux-surface
coordinate $r=0.5a$ was chosen, where $a$ is the tokamak minor radius. The
normalized density and temperature gradient profiles, for a general profile
$A(r)$, are defined as $R/L_{A}=-R\partial_{r}\textrm{ln}(A(r))$, which can be
calculated using the profile given in Eq. (2) of Ref. GorlerBenchmark . The
values of the pertinent simulation parameters are listed in Table II.

Table I. Relevant GENE normalizations and definitions GeneCode .

$\rho^{*}$ | $\rho_{i}/R$
---|---
$c_{s}$ | $\sqrt{T_{i}/m_{e}}$
$t$ | $R/c_{s}$
$\omega$ | $c_{s}/R$
$\gamma$ | $c_{s}/R$
$Q_{\textrm{gB}}$ | $c_{s}n_{e}T_{e}(\rho^{*})^{2}$
$\chi_{\textrm{gB}}$ | $\rho_{i}^{2}c_{s}/R$
$\phi$ | $e\phi/(T_{e}\rho^{*})$
$\nu_{c}$ | $\pi\textrm{ln}\Lambda e^{4}n_{e}R/(2^{3/2}T_{e}^{2})$
$\Lambda$ | $24-\textrm{ln}(\sqrt{10^{13}\cdot n_{e}}/(10^{3}T_{e}))$

Table II. GENE simulation parameters.

$R(\textrm{m})$ | $1.67$
---|---
$n_{i,e}(10^{19}\textrm{m}^{-3})$ | $4.66$
$T_{i,e}(\textrm{keV})$ | $2.14$
$B_{\phi}(\textrm{T})$ | $2.0$
$r/a$ | $0.5$
$a/R$ | $0.36$
$\rho^{*}$ | $0.001413$
$\beta$ | $1\textrm{e}{-4}$
$m_{i}/m_{p}$ | $1.0$
$m_{e}/m_{i}$ | $5.4462\textrm{e}{-4}$
$R/L_{T_{i}}$ | $0$
$R/L_{T_{e}}$ | $6.96$
$R/L_{n_{i,e}}$ | $2.22$
$q$ | $1.41$
$\hat{s}$ | $0.837$
$\nu_{ei}$ | $0.106875$

For the linear ITG case, the grid resolution was taken with $32$ grid points
in the radial dimension, $x$, and $16$ grid points in the parallel spatial
dimension, $z$. The GENE radial coordinate $x$ corresponds to the flux-surface
coordinate $r$ for the case of a circular geometry. In Fig. 1 the flux-tube
GENE benchmark result from Ref. GorlerBenchmark  is marked as “Görler” and the
corresponding “low” ($32\times 8$) and “high” ($64\times 16$) velocity grid
($v_{\parallel}\times\mu$) resolution simulations have been plotted
collectively. There is good agreement with the benchmark case in the ITG
range. The intermediate scale is well resolved in the “low v-res” case, and
the same grid resolution was used for the nonlinear simulation, but with the
radial grid resolution increased to $192$ gridpoints. The perpendicular box
size was reduced from $125\rho_{i}\times 125\rho_{i}$ at the ITG turbulence
scale to $6\rho_{i}\times 3\rho_{i}$ at the ETG turbulence scale, where the
original radial extent of the flux-tube domain was increased from $3\rho_{i}$
to $6\rho_{i}$ to allow for the full formation of the ETG mode streamers.

The electrostatic portion of the radial heat flux, $\langle Q_{ES}\rangle$,
for electrons is shown approaching a statistically steady state in time in
Fig. 2 for the collisionless, $6\rho_{i}\times 3\rho_{i}$, nonlinear, full-
spectrum, electron-scale case. The heat flux is normalized to $Q_{gB}$, the
gyroBohm normalization given in Table I, and the angled brackets,
$\langle...\rangle$, denote a flux-tube volume average. The heat flux is
determined in GENE as GeneCode ,

$Q_{ES}=\int d^{3}v\left(\frac{1}{2}mv^{2}\right)\bm{v}_{E\times
B}\cdot\hat{\bf r}\;\delta f\;,$ (1)

where $\delta f$ is the distribution function perturbation, $\bm{v}_{E\times
B}=-(\nabla\phi\times\bm{B})/B^{2}$ is the $\bm{E}\times\bm{B}$ drift, $\phi$
is the electrostatic potential perturbation, and $m$ and $v$ are the particle
mass and velocity. The normalization for $\phi$ is given in Table I. One can
see the shift from the early, nonlinear state characterized by radially-
elongated electrostatic potential streamers in Fig. 2 to the later state of
Fig. 2 where zonal flows have become dominant. It is during this phase that
intermediate-scale zonal flows grow slowly into the final quasi-saturated
state.

The initial condition for $\phi$ was realistically peaked about the most
unstable mode as shown in Fig. 3. This allowed for a transition from the
high-$k_{y}$ ETG turbulence regime to the intermediate scale where ZF
generation is expected to be stronger Chen_ETG_ZF . Such zonal flow generation
is not present in toroidal electron fluid theories LinFluidETG ; ChenFluidETG
. An inverse-cascade can clearly be seen between Figs. 3 and 3. This initial
saturation is discussed further for the single-mode simulation results
presented in Sec. III.2, and for the well-converged, collisional, full-
spectrum simulation results presented in Sec. IV. The convergence tests for
finding an optimal nonlinear box size are detailed in Appendix A.

##  III Single-Mode Analysis

###  III.1 Zonal Flow Generation Mechanism

It was shown theoretically in Ref. Chen_ETG_ZF  that intermediate-scale zonal
flow may play a role in the nonlinear saturation of ETG turbulence in tokamak
plasmas. The equations of Ref. Chen_ETG_ZF  are briefly reviewed here, and the
nonlinear single-mode simulation results are then presented in Sec. III.2 with
comparison to the theory. The theoretical model takes the standard gyrokinetic
equation FriemanChenGyro  for electrons and the quasineutrality condition in
the intermediate-scale limit $k_{\perp}^{2}\rho_{e}^{2}\ll 1\ll
k_{\perp}^{2}\rho_{i}^{2}$. For long-wavelength ETG modes, one can generally
assume that the growth rate is much smaller than the real frequency,
$|\gamma/\omega|\ll 1$, which allows for a weak-turbulence analysis. These
modes also satisfy the relation $|\omega_{d}|\ll|\omega|,|\omega_{t}|$, for
$\omega_{t}$ and $\omega_{d}$ the transit and magnetic drift frequencies
respectively. Local approximations KimHortLocal  are then assumed for
$\omega_{d}$ and $\omega_{t}$ in formulating a kinetic electron model.

The evolution of a single ETG mode can then be derived in terms of a nonlinear
Schrodinger equation (NLSE) Chen_ETG_ZF ,

$\displaystyle[i(\partial_{t}-\gamma_{n})-b_{n}k_{\theta}^{2}\rho_{e}^{2}\hat{s}^{2}\theta_{k}^{2}-\frac{c_{n}}{k_{\theta}^{2}\rho_{e}^{2}\hat{s}^{2}}\frac{\partial^{2}}{\partial\theta_{k}^{2}}]A_{n}(\theta_{k})$
(2) $\displaystyle=\frac{-ik_{\theta}\rho_{e}\hat{s}}{\sqrt{2\pi}}\int
d\vartheta_{k}\vartheta_{k}A_{z}(\vartheta_{k})A_{n}(\theta_{k}-\vartheta_{k}),$

with $A_{z}$ and $A_{n}$ the amplitudes of the ZF and ETG modes, $\gamma_{n}$
the ETG mode growth rate, $\theta_{k}$ and $\vartheta_{k}$ the tilting angle
(defined by $\theta_{k}=k_{x}/(k_{y}\hat{s})$) in the flux-tube limit
BeerThesis  for $k_{y}$ of the ETG mode), $k_{\theta}$ the poloidal
wavenumber, and the terms including $b_{n}$ and $c_{n}$ coming from the
frequency mismatch and plasma nonuniformity corrections respectively. $b_{n}$
and $c_{n}$ are specifically associated with the linear ETG dynamics as
explained in Ref. Chen_ETG_ZF . The nonlinear term under the integral in Eq.
(2) describes a Navier-Stokes type coupling due to E$\times$B shearing
effects. This coupling is $\mathcal{O}((k_{\theta}\rho_{e})^{-2})$ stronger
than the usual Hasegawa-Mima type coupling in the fluid limit LinFluidETG ;
ChenFluidETG . This stronger coupling results in a stronger regulation of ETG
turbulence by zonal flows and also leads to a reduced threshold for ZF
excitation by intermediate-scale ETG modes. The threshold condition is
described further in Eq. (4).

The description of zonal flow dynamics is given by the equation Chen_ETG_ZF ,

$\displaystyle[\partial_{t}+\gamma_{z}(1+d_{z}k_{\theta}^{2}\rho_{e}^{2}s^{2}\theta_{k}^{2})]\chi_{z}A_{z}(\theta_{k})$
(3) $\displaystyle=\sqrt{\frac{\pi}{2}}(k_{\theta}\rho_{e}s)\theta_{k}\int
d\vartheta_{k}\vartheta_{k}^{2}[A_{n}(\vartheta_{k})A_{n}^{*}(\vartheta_{k}-\theta_{k})a_{n}^{*}$
$\displaystyle-A_{n}(\vartheta_{k}+\theta_{k})A_{n}^{*}(\vartheta_{k})a_{n}].$

The nonlinear term under the integral in Eq. (3) is related to the Reynolds
stress of the ETG modes. The ZF damping rate is given by $\gamma_{z}\simeq
3\nu_{ee}/(|\omega_{*e}|\sqrt{\epsilon})$, with $\nu_{ee}$ the electron-
electron collision frequency and $\omega_{*e}$ the electron diamagnetic drift
frequency. The total electric susceptibility is defined as
$\chi_{z}=\tau+(1+1.6q^{2}/\sqrt{\epsilon})k_{\theta}^{2}\rho_{e}^{2}s^{2}\theta_{k}^{2}/2$
with $\tau=T_{e}/T_{i}$. The $d_{z}$ term in Eq. (3) represents a
gyrodiffusive correction which helps suppress short wavelength zonal flows
GyroDiffDampZF .

The $a_{n}$ term describes the parallel correlation of the ETG turbulence and
is defined as $a_{n}(\theta_{k},\vartheta_{k})=\int
d\eta\langle\tilde{\Phi}_{n}^{*}(\eta,\vartheta_{k})v_{\perp}^{2}\delta\tilde{H}_{n}(\eta,\vartheta_{k}+\theta_{k})\rangle_{v}$,
where $\eta$ is the extended poloidal angle. Here, $\delta\tilde{H}_{n}$, is
the non-adiabatic part of the distribution function perturbation. Keeping
$a_{n}$ non-local in tilting angle in order to take into account ballooning
effects leads to a parallel decoupling of ETG modes and therefore ZFs as well,
where $a_{n}=a_{n}(0,0)\textrm{exp}(-\theta_{k}^{2}/2\bar{\eta}^{2})$ for
$\bar{\eta}=(\int d\eta\tilde{\Phi}^{*}\eta^{2}\tilde{\Phi})^{1/2}$ the
parallel-mode-averaged potential Chen_ETG_ZF .

Eqs. (2) and (3) are taken together as the NLSE model. The numerical solution
of the NLSE model with a single ETG mode and a range of ZF modes gives an
evolution of ETG and ZF modes that can be described by three specific stages
Chen_ETG_ZF . The initial stage involves uninhibited exponential growth of the
ETG mode to a threshold point at which the radial beating of the ETG drives ZF
growth as described by equation (3). As the ZF modes grow, they lead to radial
dispersion of the initial ETG wave packet and the creation of sidebands. These
sidebands then drive more zonal flows via a modulational instability in the
second stage. Once the zonal flow grows to appreciable levels in comparison
with the ETG mode, the nonlinear interaction in equation (2) acts to saturate
the ETG mode. In the final stage, the linear growth rate of the ETG mode
becomes negligible and the NLSE model then results in slow, algebraic growth
for the zonal flow Chen_ETG_ZF . This slow growth has been observed in
previous gyrokinetic electron-scale simulations ParkerZF ; ColyerETG .

A threshold condition for the ZF excitation can be calculated analytically
Chen_ETG_ZF  by considering a simple four-wave model for a single zonal flow
mode, an ETG pump mode, and two ETG sideband modes. Narrow-band, rectangular
functions are used to describe the ZF and ETG modes,
$A_{z}\Pi[(\theta_{k}-\theta_{z})/W]$ and
$A_{0}\Pi(\theta_{k}/W)+A_{+}\Pi[(\theta_{k}-\theta_{z})/W]+A_{-}\Pi[(\theta_{k}+\theta_{z})/W]$
respectively. Here $W$ is the full-width of the modes, $\theta_{z}$ is the ZF
wavenumber in terms of tilting angle, and $A_{z},A_{0},$ and $A_{\pm}$ are the
ZF, ETG pump, and ETG sideband mode amplitudes. Substituting these functions
into Eqs. (2) and (3) with the assumption of no plasma nonuniformities
($c_{n}=0$) and a steady state pump amplitude for simplicity, one obtains the
critical threshold condition Chen_ETG_ZF ,

$W^{2}|A_{0,c}^{2}|=\frac{(\Delta^{2}+\gamma_{s}^{2})\gamma_{z}(1+d_{z}(k_{\theta}\rho_{e}\hat{s}\theta_{z})^{2})\chi_{z}}{(k_{\theta}\rho_{e}\hat{s}\theta_{z})^{4}[\gamma_{s}\textrm{Re}(a_{n})+\Delta\textrm{Im}(a_{n})]}.$
(4)

Here, $\Delta=\textrm{Re}(b_{n})(k_{\theta}\hat{s}\theta_{z}\rho_{e})^{2}$
represents the frequency mismatch of the pump and sidebands, and
$\gamma_{s}=\gamma_{n}+\textrm{Im}(b_{n})(k_{\theta}\hat{s}\theta_{z}\rho_{e})^{2}$
is the growth rate of the sidebands. This threshold condition for ZF
excitation by intermediate-scale ETG modes is
$\mathcal{O}(k_{\theta}^{2}\rho_{e}^{2})$ lower than the condition found in
the fluid approximation ChenFluidETG , which would then lead to more effective
ZF generation at intermediate scales.

Figure 4: Plots of collisionless single ETG mode evolution for (a)
$k_{y}\rho_{i}=6.36$ and (b) $k_{y}\rho_{i}=30$ with respective growth rates
$\gamma\approx 7.037$ and $7.015$. The strongest four ZF modes in the late
time are plotted for (c) the $k_{y}\rho_{i}=6.36$ case, while (d) shows the
excitation of a larger range of ZF modes in the $k_{y}\rho_{i}=30$ case.
Figure 5: Plot of the sum of zonal potential amplitudes at final times,
$\sum|\phi_{z}|=\sum_{k_{x}}|\phi_{k_{x},0}|$, as driven by a single ETG mode.
The fourth and seventh points correspond to the ETG modes from Fig. 4.

Figure 6: Plot of total (a) ETG and (b) ZF mode potential as an integral over
$k_{x}$ in the collisional, Gaussian-$k_{x}$ ETG case for
$k_{y}\rho_{i}=6.36$. Markers have been added to match the spectral snapshots
of Fig. 7.

Figure 7: Radial-spectrum snapshots of (a) the single ETG mode and (b) ZF as a
function of tilting angle for the $k_{y}\rho_{i}=6.36$ case. Snapshot times
match the markers in Fig. 6.

The ZF modes are initially excited in a range of $\theta_{z}$ values. As the
system evolves to the quasi-saturated state and the ETG mode is suppressed,
the ZF spectrum narrows towards the most easily driven mode. Comparing the
exponential $\theta_{k}$ dependence of the $a_{n}$ term to the algebraic form
of the $d_{z}$ term, one finds that the parallel decoupling is largely
responsible for minimizing the intermediate-scale threshold condition at low
$\theta_{z}$. If one considers the temporal evolution of the four-wave model,
a fixed-point solution can be found with constant $A_{z}$ and sufficiently low
ZF damping rate, where the ZF and ETG mode amplitudes are given as Chen_ETG_ZF
,

$W^{2}|A_{z,p}|^{2}=\frac{\pi(\delta^{2}-\Delta\delta-\gamma_{n}\gamma_{s})}{(k_{\theta}\hat{s}\theta_{z}\rho_{e})^{2}},$
(5)

and,

$W^{2}|A_{0,p}|^{2}=\frac{\chi_{z}\gamma_{z}(1+d_{z}k_{\theta}^{2}\hat{s}^{2}\theta_{z}^{2}\rho_{e}^{2})[\gamma_{s}^{2}+(\delta-\Delta)^{2}]}{(k_{\theta}\hat{s}\theta_{z}\rho_{e})^{4}[(\Delta-\delta)\textrm{Im}(a_{n})+\gamma_{s}\textrm{Re}(a_{n})]}.$
(6)

Here, $\delta=\Delta\gamma_{n}/(\gamma_{s}+\gamma_{n})$ represents the
amplitude oscillation frequency of ETG modes due to the nonlinear ETG-ZF
coupling. The single ETG mode spectrum then continues to fluctuate in $k_{x}$
while the ZF mode reaches a constant, steady state Chen_ETG_ZF . Additionally,
$|A_{0,p}|^{2}$ is then proportional to $\gamma_{z}$ in the saturated state,
and therefore the collision frequency, while $|A_{z,p}|^{2}$ is not. While
these saturation estimates are only valid for a single ZF mode, as the ETG
turbulence saturates and the ZF spectrum narrows due to the threshold
condition, $\theta_{z}$ of the most optimally-driven mode can be used to
estimate the ETG saturation level.

###  III.2 Single-Mode Simulation Results

We now compare the NLSE model (Eqs. (2) and (3)) to gyrokinetic simulation
results. Collisionless, nonlinear, ETG simulations were carried out where a
single unstable ETG mode ($k_{x}=0$) and all zonal flow modes ($k_{y}=0$) are
retained. This fairly accurately describes the dynamics of the NLSE model. All
results presented in this section are averaged over $z$. The ETG growth rate
spectrum with respect to $k_{y}\rho_{i}$ can be utilized here to illustrate
the NLSE model dynamics. As seen in Fig. 1, the ETG growth rate spectrum is
symmetric around the most unstable mode, so one can choose to compare the
evolution of a pair of ETG modes with similar growth rates, where one mode has
a $k_{y}\rho_{i}$ value in the intermediate-scale range and the other mode has
a larger $k_{y}\rho_{i}$ value outside of that range. Then the ZF drive of the
two modes can be compared to verify the expectations from the NLSE model.

The $k_{y}\rho_{i}=6.36$ mode with a growth rate of $\gamma\approx 7.037$ and
the $k_{y}\rho_{i}=30$ mode with a similar growth rate of $\gamma\approx
7.015$ are taken here for comparison. Fig. 4(a) shows the time evolution of
the $k_{y}\rho_{i}=6.36$ mode, while Fig. 4(b) shows the time evolution of the
$k_{y}\rho_{i}=30$ mode. Fig. 4(c) shows the time evolution of the four
strongest zonal flow modes at the final time step for the intermediate-scale
case, whereas Fig. 4(d) shows a large range of ZF modes excited in the
high-$k_{y}$ case in order to illustrate a difference in the zonal flow
response between the two cases. One can see that initially both ETG modes grow
exponentially at similar rates until a threshold is reached, at which point
zonal flows are excited. For the intermediate-scale ETG mode, this phase is
followed by an algebraically-growing long wavelength ZF phase in which the ZF
modes gradually reach a steady state value. The high-$k_{y}$ results show no
slowly-growing ZF phase at late times. This difference in ZF generation in the
late stage is consistent with the threshold condition given in Eq. (4). The
intermediate-scale ETG mode continues to slowly drive zonal flows as it is
suppressed to lower levels, whereas the high-$k_{y}$ ETG mode does not.

The peak level of the ETG mode is much lower for the high-$k_{y}$ case than
for the intermediate-scale case. This result is not expected from the NLSE
model as the shearing of the ETG mode by the wave-wave coupling should be
stronger in the intermediate scale than at higher $k_{y}$. Zonal flows are
also generated earlier in the high-$k_{y}$ case, indicating a lower threshold
initially, which is inconsistent with the NLSE model. One noticeable
difference between the single-mode GENE simulations and the NLSE model is that
the NLSE model only includes the ZF shearing suppression mechanism, whereas
the single-mode GENE simulations include other saturation mechanisms.
Comparing the zonal response between the two cases, it is found that the
initial ZF excitation shown in Fig. 4(d) is much more abrupt, possibly
indicative of a secondary instability DorlandJenkoETG1 ; DorlandJenkoETG2 ;
ETGFourSec . The change in ETG-ZF dynamics in the single-mode results is found
to occur near $k_{y}\rho_{i}\sim 15$. This difference in behavior likely
indicates the reason for the transition to the intermediate scale mentioned
previously in Sec. II, and this is further discussed in comparison to the
full-spectrum simulation results presented in Sec. IV.

Fig. 5 shows the sum of all ZF mode amplitudes, $\sum|\phi_{z}|$, as a
function of $k_{y}\rho_{i}$. Each value of $k_{y}\rho_{i}$ represents
initializing with a different unstable ETG mode. The sum is taken at the final
simulation time, where the ZF mode amplitudes are nearing steady state levels.
The notable region of ZF generation is clearly seen to be in the intermediate-
scale range, as expected by the NLSE model. Shorter ETG mode wavelengths
correspond to weaker zonal flows at late times, in agreement with fluid ETG
models. In addition, the drop-off at long ETG mode wavelengths is reasonable
due to trapped electron effects at this wavenumber range ChenCTEM3 . A
validation of this expectation for the full-spectrum simulations is provided
in Appendix A.

The unstable ETG mode is shown to be suppressed at late times in Figs. 4(a)
and 4(b). The total amplitude, $\langle|\phi|\rangle=(\int
d\theta_{k}|\phi_{k}|^{2})^{1/2}$, of the ETG mode is small in comparison to
the zonal flow amplitude. In contrast, the NLSE model simulation results given
in Fig. 3 of Ref. Chen_ETG_ZF  show that the total ETG and ZF mode amplitudes,
$\langle|\phi|\rangle$, are of similar strength and fluctuating in the late
stage. One reason to expect the strong ETG suppression in the gyrokinetic
simulations is the lack of collisionality which would damp the zonal flow due
to the collisional dependence of the $\gamma_{z}$ term in Eq. (3).
Additionally, the NLSE model assumes a Gaussian radial spectrum for the ETG
mode, while the single-mode flux-tube GENE simulations take $k_{x}=0$
initially. The globally-Gaussian radial distribution of the ballooning modes
would lead to more radial ETG mode overlap, which would then drive more ZF
generation as predicted by Eq. (3).

Figure 8: Plot in time of the average dimensionless radial wavenumber for the
total ETG and ZF spectra in the quasi-saturated stage.

We were able to obtain more physical results for the $k_{y}\rho_{i}=6.36$ case
in which the ETG and ZF modes fluctuate about similar steady state values due
to their nonlinear coupling. These results were achieved by including the
physical collisionality taken in Sec. IV, such that $\gamma_{z}=0.014$, and by
initializing the ETG mode with a Gaussian $k_{x}$ and $z$ spectrum such that
$\phi(k_{x},z,t=0)\sim e^{-(k_{x}^{2}+z^{2})/8\pi^{2}}$. The corresponding
$\langle|\phi|\rangle$ for the ETG and ZF modes are plotted in Fig. 6. One can
see the strong drive of the zonal flow, as well as the late-stage fluctuations
of both ETG and ZF modes. Fig. 7 shows that over time the Gaussian radial
spectrum of the ETG mode is broadening into sidebands, while the radial
spectrum of the zonal flow modes narrows from a broad distribution to a peak
at a final, low-$\theta_{z}$ mode number, as predicted by the NLSE model
Chen_ETG_ZF . These results suggest that one should perform global simulations
to see results most consistent with the NLSE model.

The improved single-mode results for the $k_{y}\rho_{i}=6.36$ case showed
large $\theta_{k}$-averaged fluctuations for the ETG mode compared to more
fixed ZF fluctuations. These results agree with the expectations of Eqs. (5)
and (6). The fluctuations can be seen in Fig. 8, where the average
dimenionsless radial wavenumber, $\langle|\theta_{k}|\rangle=\left(\int
d\theta_{k}\theta_{k}^{2}|\phi_{k}|^{2}\right)^{1/2}/\langle|\phi_{k}|\rangle$,
is plotted for the ETG and ZF modes respectively. The ratio of the total
absolute amplitude of ZF to ETG modes is given in Fig. 9 as a function of
$\gamma_{z}$ and is consistent with the trend from the NLSE model. The late-
stage behavior of the ETG and ZF $k_{x}$-spectra, as shown in Figs. 7 and 8,
and the collisional behavior of the mode amplitude ratio agree well with the
late-time behavior reported in the electron-scale MAST simulations of Ref.
ColyerETG .

Figure 9: Ratio of total absolute amplitudes of ZF to ETG modes with varying
collisionality. Results are taken at the final timestep. The reference value
corresponds to the fourth point, $\gamma_{z}=0.014$.

##  IV Full-spectrum Simulation Results

The full-spectrum nonlinear simulation results are presented here and the
intermediate-scale zonal flow generation mechanism is further investigated.
Including multiple toroidal modes results in a final quasi-saturated heat flux
characterized by richer turbulent interactions. Fig. 10(a) shows the time
history of the heat flux for the well-converged $24\rho_{i}\times 3\rho_{i}$
case with collisionality. The four strongest zonal flow modes at the final
time are presented in Fig. 10(b). The value of the normalized electron-ion
collision frequency used is $\nu_{ei}=0.106875$. This frequency is defined as
$\nu_{ei}=4v_{te}\nu_{c}/R$, where $v_{te}=\sqrt{T_{e}/m_{e}}$ and $\nu_{c}$
is the collision frequency given in Table I. The self-adjoint form of the
standard Landau-Boltzmann collision operator is used. Realistic collisionality
allows for ZF damping when reaching a final state, and the simulation was
carried out to a sufficiently long non-dimensional time, $t/(R/c_{s})=90$, to
ensure that a quasi-saturated steady state in $\langle Q_{ES}\rangle$ is
achieved. The convergence with respect to box size is discussed in Appendix A.

Figure 10: Full-spectrum collisional nonlinear simulation results showing (a)
the electron heat flux time-evolution and (b) the time evolution of the four
strongest final ZF modes. Data is averaged over $z$. Figure 11: Time trace of
full-spectrum ETG modes ranging from the longest mode to the most unstable
mode. Data is averaged over $z$.

Fig. 11 shows the evolution of various ETG modes ranging from the longest
wavelength to the most unstable mode. The shorter wavelength modes saturate
very quickly to negligible levels, in agreement with the single-mode results,
and are omitted. One can observe that the intermediate-scale modes saturate
the slowest and reach the highest levels. During the period of intermediate-
scale ETG mode growth, the ZF modes shown in Fig. 10(b) are driven
exponentially by the radial beating of ETG modes, as well as by the
modulational instability. Once the ETG modes reach a quasi-saturated state,
the ZF modes continue to grow slowly in agreement with the single-mode
simulation results. Considering the findings of the NLSE model, the single-
mode simulations results of Sec. III.2, and the full-spectrum simulation
results discussed here, it is the intermediate-scale ETG modes which are most
responsible for driving ZF mode growth into the late stage.

Figure 12: Comparison of quasilinear and nonlinear heat flux spectra for the
$24\rho_{i}\times 3\rho_{i}$ full-spectrum electron-scale case. NL and QL
stand for nonlinear and quasilinear respectively.

The heat flux spectrum for the full-spectrum case is shown in Fig. 12,
alongside a quasilinear saturation estimate. The quasilinear estimate of the
heat flux for a single $k_{y}$ mode is calculated as FableQuasiLin ;
LapQuasiLin ,

$Q_{k_{y}}^{QL}=A_{0}\frac{(\gamma/\langle
k_{\perp}^{2}\rangle)^{2}}{|\phi_{0,k_{y}}(0)|^{2}}Q_{k_{y}}^{\textrm{lin.}},$
(7)

with $Q_{k_{y}}^{\textrm{lin.}}$ representing the linear simulation results
for the heat flux and $\phi_{0,k_{y}}(0)$ the linear electrostatic potential
at $k_{x}=z=0$. $A_{0}$ represents a constant of proportionality, and $\langle
k_{\perp}^{2}\rangle$ is the ballooning-angle-averaged perpendicular
wavenumber, defined as LapQuasiLin ,

$\langle
k_{\perp}^{2}\rangle=\frac{\sum_{k_{x}}\int(g^{xx}k_{x}^{2}+2g^{xy}k_{x}k_{y}+g^{yy}k_{y}^{2})|\phi_{k}(z)|^{2}dz}{\sum_{k_{x}}\int|\phi_{k}(z)|^{2}dz}.$
(8)

Here, $\phi_{k}$ represents a Fourier mode of the electrostatic potential
perturbation, and $g^{\mu\nu}=\nabla\mu\cdot\nabla\nu$ gives the metric tensor
coefficients of the GENE coordinate system. A sum over all $k_{x}$ values was
required for the drop-off of saturation levels at low $k_{y}\rho_{i}$.

This model of mode saturation describes a balance of the unstable growth of
the instability with turbulent diffusion based on a mixing-length estimate
JenkoSat1 ; JenkoSat2 ; BourdSat , and gives notable agreement with the
nonlinear heat flux spectrum in the intermediate scale. The disparity between
the quasilinear model and the nonlinear heat flux spectrum is greatest for ETG
modes with wavenumbers higher than $k_{y}\rho_{i}\approx 12$. This suggests
other, stronger saturation mechanisms for these modes. While new effects, such
as toroidal inverse-cascading LinFluidETG ; ChenFluidETG ; KimETG , may play a
role in the full-spectrum case, the difference in the heat flux spectra also
agrees with the transition in ETG-ZF dynamics found to occur around
$k_{y}\rho_{i}=15$ for the single-mode results. The more abrupt ZF response
and quicker saturation of the higher $k_{y}$ modes may be consistent with
secondary instability theory DorlandJenkoETG1 ; DorlandJenkoETG2 , where a
$|\phi_{n}|\sim\gamma/\langle k_{\perp}^{4}\rangle$ saturation model for ETG
mode amplitudes predicts a steeper drop-off near $k_{y}\rho_{i}\approx 12$
than that of the quasilinear mixing-length estimate.

Finally, we present a comparison of turbulent and neoclassical transport
levels at both ion and electron scales. Because the electron-scale case takes
the ion temperature gradient to zero, we can compare the electron-scale
thermal diffusivity to that of the ion-scale ITG case with adiabatic electrons
to understand the importance of regulation by zonal flows at each scale. In
units normalized to the specific species of interest, the thermal diffusion
coefficients due to electrostatic turbulence are
$\langle\chi_{\textrm{ES}}\rangle_{i}=0.7\rho_{i}^{2}v_{Ti}/L_{Ti}$ for the
ion-scale ITG case and
$\langle\chi_{\textrm{ES}}\rangle_{e}=2.8\rho_{e}^{2}v_{Te}/L_{Te}$ for the
electron-scale ETG case, where $v_{Ts}$ is the thermal velocity,
$\sqrt{2T_{s}/m_{s}}$, for a species $s$. This suggests that the ETG-driven
zonal flows don’t regulate ETG turbulence as strongly as the isomorphic
counterpart ITG turbulence is regulated by ITG-driven zonal flows.

The neoclassical transport values were calculated using GENE for both the ion
and electron scale cases. Given in units of $\chi_{gB}$ from Table I, the
neoclassical thermal diffusivites are
$\langle\chi_{\textrm{neo}}\rangle_{i}=0.14\chi_{\textrm{gB}}$ and
$\langle\chi_{\textrm{neo}}\rangle_{e}=0.004\chi_{\textrm{gB}}$. The
neoclassical values are in close agreement with the theoretical expectation
that $\chi_{i}=\sqrt{m_{i}/m_{e}}\chi_{e}$, and are negligible compared to the
turbulent thermal diffusivities,
$\langle\chi_{\textrm{ES}}\rangle_{i}=6.95\chi_{\textrm{gB}}$ and
$\langle\chi_{\textrm{ES}}\rangle_{e}=0.328\chi_{\textrm{gB}}$. The late-time
heat flux spectrum peaks in the intermediate scale at $k_{y}\rho_{i}=10.6$
with $\langle Q_{ES}\rangle_{e}=0.66Q_{\textrm{gB}}$ and drops off to $\langle
Q_{ES}\rangle_{e}=0.10Q_{gB}$ and $0.11Q_{gB}$ for $k_{y}\rho_{i}=4.24$ and
$16.96$ respectively. These values are in good agreement with the theoretical
expectation that $Q/Q_{gB}\sim\mathcal{O}(0.01)-\mathcal{O}(0.1)$ (in the
units of Table I) for the intermediate-scale ETG modes Chen_ETG_ZF .

##  V Discussion

We have shown, using the single-mode nonlinear simulations, that the NLSE
model Chen_ETG_ZF  accurately describes the zonal flow generation mechanism by
intermediate-scale ETG modes and that it provides a theoretical understanding
for the slow growth of long-wavelength zonal flows into the long-term quasi-
saturated state. As the NLSE model considers only a single ETG mode for a
practicable analysis, one cannot say conclusively that the same is true of the
full-spectrum nonlinear results. However, in the full-spectrum case the
high-$k_{y}$ ETG modes are quickly saturated by a stronger ZF response as
compared to the intermediate-scale ETG modes. The intermediate-scale ETG modes
then drive exponential ZF mode growth initially, and slow, algebraic ZF mode
growth as they are suppressed in the late stage. This result is in good
agreement with the NLSE model for intermediate-scale ETG-ZF dynamics, as well
as various long time, saturated electron-scale ParkerZF ; ColyerETG  and
multiscale HowardLmode ; HollandHmode  flux-tube simulations.

Figure 13: Convergence of the electrostatic electron heat flux for the various
box size cases compared to the original $6\rho_{i}\times 3\rho_{i}$ case.

Figure 14: (a) ETG and (b) four strongest final ZF mode time traces for the
$12\rho_{i}\times 6\rho_{i}$, large-y box size. All data is averaged over $z$.

The final transport levels for the full-spectrum case are in similar ranges
found in thorough electron-scale CBC benchmarks which compare well with
experimental observations Nevins06 . As the zonal flows are driven at long
electron-scale wavelengths, multiscale effects could become important and ETG-
driven zonal flows may have an effect on ion-scale turbulence. This effect
where intermediate-scale zonal flows contribute to ion-scale turbulence
suppression has been reported in large multiscale simulations HowardLmode ;
HollandHmode .

## Acknowledgements

This research was supported by the SciDAC-4 project High-fidelity Boundary
Plasma Simulation funded by the U.S. Department of Energy (U.S. DOE) Office of
Science under grant DE-SC-000801. Research was also carried out under U.S. DOE
grant DE-FG02-08ER54954.

##  Appendix A Nonlinear Convergence Tests

This section details a “full-spectrum” nonlinear simulation box-size
convergence study. The perpendicular box size is varied in terms of the basic
$6\rho_{i}\times 3\rho_{i}$ electron-scale box size shown in Fig. 2. The
collisionality was set to the reference value discussed in Sec. IV. Four
perpendicular domain sizes are presented: $6\rho_{i}\times 3\rho_{i}$,
$12\rho_{i}\times 3\rho_{i}$, $24\rho_{i}\times 3\rho_{i}$, and
$12\rho_{i}\times 6\rho_{i}$. These cases consider the importance of correctly
resolving the ETG streamer lengths and the longest wavelength zonal flow
modes. Additionally, the inclusion of longer ETG mode wavelengths is
considered in the $12\rho_{i}\times 6\rho_{i}$ case to verify the findings of
Section III.2. The number of radial gridpoints was increased in each
simulation to retain the original resolution.

The electron heat flux for each case is shown over time in Fig. 13. One can
clearly see that the increase in radial dimension is necessary to correctly
resolve the heat flux. The late-stage zonal flows of the $6\rho_{i}\times
3\rho_{i}$ case, as shown in Fig. 3(c), are peaking at the longest mode
allowable and the box size must be increased to correctly resolve the longest
modes. Allowing for longer wavelength ZF modes leads to stronger regulation of
the heat flux as seen in Fig. 13.

The time evolution of the ETG modes for the $12\rho_{i}\times 6\rho_{i}$ case
is shown in Fig. 14(a), and the time evolution of the four strongest ZF modes
at the final time is shown in Fig. 14(b). These results are qualitatively
similar to the $24\rho_{i}\times 3\rho_{i}$ case shown in Figs. 10(b) and 11.
In this new case, the longest wavelength ETG mode, $k_{y}\rho_{i}=1.06$, grows
to the highest level. However, it can be seen in Fig. 14(b) that from non-
dimensional times 5-10 $t/(R/c_{s})$, when the longest ETG mode is dominant,
the zonal flows are already in the final, slowly growing stage. This result
indicates that the strongest ZF modes are largely being affected by the
intermediate-scale ETG modes, not the longest wavelength ETG mode, and
confirms the results found in Sec. III.2 which showed little zonal flow
generation outside the intermediate-scale range. As increasing $L_{y}$ from
the original $3\rho_{i}$ size had no effect on the final quasi-saturated
state, the largest $L_{x}$ case considered, $24\rho_{i}\times 3\rho_{i}$, was
chosen for the full-spectrum investigation discussed in Sec. IV.

## References

  * [1] Chen H. et al 2021 Nucl. Fusion 61 066017
  * [2] Batchelor D. A. et al 2007 Plasma Sci. Technol. 9 312
  * [3] ITER Physics Expert Group on Confinement and Transport et al 1999 Nucl. Fusion 39 2175
  * [4] Horton W. 1999 Rev. Mod. Phys. 71 735
  * [5] Doyle E.J. et al 2007 Nucl. Fusion 47 S18
  * [6] Lee Y.C., Dong J.Q., Guzdar P.N., and Liu C.S. 1987 Phys. Fluids 30 1331
  * [7] Jenko F., Dorland W., Kotschenreuther M., and Rogers B.N. 2000 Phys. Plasmas 7 1904
  * [8] Dorland W., Jenko F., Kotschenreuther M., and Rogers B.N. 2000 Phys. Rev. Lett. 85 5579
  * [9] Jenko F. and Dorland W. 2002 Phys. Rev. Lett. 89 225001
  * [10] Nevins W.M., Candy J., and Cowley S. et. al. 2006 Phys. Plasmas 13 122306
  * [11] Ren Y. et al 2017 Nucl. Fusion 57 072002
  * [12] Gierson B.A., Stabler G.M., Solomon W.M. et al 2018 Phys. Plasmas 25 022509
  * [13] Ryter F. et al 2019 Nucl. Fusion 59 096052
  * [14] Kiefer C.K. et al 2021 Nucl. Fusion 61 066035
  * [15] Howard N.T. et al 2016 Phys. Plasmas 23 056109
  * [16] Holland C. et al 2017 Nucl. Fusion 57 066043
  * [17] Parker S.E. et al 2006 AIP Conf. Proc 871 193
  * [18] Colyer G.J. et al 2017 Plasma Phys. Control. Fusion 59 055002
  * [19] Hasegawa A. and Mima K. 1977 Phys. Rev. Lett. 39 205
  * [20] Lin Z., Chen L. and Zonca F. 2005 Phys. Plasmas 12056125
  * [21] Chen L., Zonca F. and Lin Z. 2005 Plasma Phys. Control Fusion 47 B71
  * [22] Kim E.J., Holland C. and Diamond P.H. 2003 Phys. Rev. Lett. 91 075003
  * [23] Dannert T. and Jenko F. 2005 Phys. Plasmas 12 072309
  * [24] Jenko F. et al 2005 Plasma Phys. Control. Fusion 47 B195
  * [25] Bourdelle C., Garbet X., Imbeaux F. et al 2007 Phys. Plasmas 14 112501
  * [26] Parker S.E. et al 1994 Phys. Plasmas 1 1461
  * [27] Diamond P.H. et al 2005 Plasma Phys. Control. Fusion 47 R35
  * [28] Jenko F. and The GENE development team, The GENE code https://genecode.org
  * [29] Greenfield C.M., DeBoo J.C., Osborne T.H., Perkins F.W., Rosenbluth M.N., and Boucher D. 1997 Nucl. Fusion 37 1215–1228
  * [30] Dimits A.M., Bateman G., Beer M.A., Cohen B.I., Dorland W., Hammett G.W., Kim C., Kinsey J.E., Kotschenreuther M., Kritz A.H., Lao L.L., Mandrekas J., Nevins W.M., Parker S.E., Redd A.J., Shumaker D.E., Sydora R., and Weiland J. 2000 Phys. Plasmas 7 969–983
  * [31] Lapillone X., Brunner S., Danner T., Jolliet S., Marinoni A., Villard L., Görler T., Jenko F., and Merz F. 2009 Phys. Plasmas 16 032308
  * [32] Görler T., Tronko N., Hornsby W.A., Bottino A., Kleiber R., Norscini C., Grandgirard V., Jenko F., and Sonnendrücker E. 2016 Phys. Plasmas 23 072503
  * [33] Rewoldt G., Lin Z., and Idomura Y. 2007 Comput. Phys. Commun. 177 775–780
  * [34] Adam J.C., Tang W.M., and Rutherford P.H. 1976 The Physics of Fluids 19 561
  * [35] Chen H.-T. and Chen L. 2018 Plasma Phys. Controlled Fusion 60 055011
  * [36] Dannert T. and Jenko F. 2005 Phys. Plasmas 12 072309
  * [37] Xiao Y. and Lin Z. 2009 Phys. Rev. Lett. 103 085004
  * [38] Chen H. and Chen L. 2019 Nucl. Fusion 59 074003
  * [39] Frieman R.A. and Chen L. 1982 Phys. Fluids 25 502
  * [40] Kim J.Y. and Horton W. 1991 Phys. Fluids B 3 1167
  * [41] Beer M.A., Ph.D., Princeton 1995 Gyrofluid Models of Turbulent Transport in Tokamaks
  * [42] Ricci P., Rogers B.N. and Dorland W. 2010 Phys. Plasmas 17 072103
  * [43] Strintzi D. and Jenko F. 2007 Physics of Plasmas 14 042305
  * [44] Chen H. and Chen L. 2022 Phys. Rev. Lett. 128 025003
  * [45] Fable E. et al 2010 Plasma Phys. Control. Fusion 52 015007
  * [46] Lapillone X. et al 2011 Plasma Phys. Control. Fusion 53 054011