# Full-F Turbulent Simulation in a Linear Plasma Device using a Gyro-Moment
Approach

B. J. Frei<EMAIL_ADDRESS>Ecole Polytechnique Fédérale de Lausanne
(EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland  Max-Planck-
Institut für Plasmaphysik, D-85748 Garching, Germany  J. Mencke  Ecole
Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015
Lausanne, Switzerland  P. Ricci  Ecole Polytechnique Fédérale de Lausanne
(EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland

###### Abstract

Simulations of plasma turbulence in a linear plasma device configuration are
presented. These simulations are based on a simplified version of the
gyrokinetic (GK) model proposed by B. J. Frei et al. [J. Plasma Phys. 86,
905860205 (2020)] where the full-F distribution function is expanded on a
velocity-space polynomial basis allowing us to reduce its evolution to the
solution of an arbitrary number of fluid-like equations for the expansion
coefficients, denoted as the gyro-moments (GM). By focusing on the
electrostatic and neglecting finite Larmor radius effects, a full-F GM
hierarchy equation is derived to evolve the ion dynamics, which includes a
nonlinear Dougherty collision operator, localized sources, and Bohm sheath
boundary conditions. An electron fluid Braginskii model is used to evolve the
electron dynamics, coupled to the full-F ion GM hierarchy equation via a
vorticity equation where the Boussinesq approximation is used. A set of full-F
turbulent simulations are then performed using the parameters of the LArge
Plasma Device (LAPD) experiments with different numbers of ion GMs and
different values of collisionality. The ion distribution function is analyzed
illustrating the convergence properties of the GM approach. In particular, we
show that higher-order GMs are damped by collisions in the high-collisional
regime relevant to LAPD experiments. The GM results are then compared with
those from two-fluid Braginskii simulations, finding qualitative agreement in
the time-averaged profiles and statistical turbulent properties.

††preprint: AIP/123-QED

##  I Introduction

Despite recent progress in the development of gyrokinetic (GK) codes, such as
`COGENT` Dorf and Dorr (2020), `Gkeyll` Hakim  _et al._ (2020), `GENE-X`
Michels  _et al._ (2022); Ulbl  _et al._ (2023) and `XGC` Chang  _et al._
(2017), extending the GK model from the core to the boundary remains
challenging since it requires dealing with a wide range of collisionality,
order-one fluctuations across various scales, complex magnetic field geometry,
steep pressure gradients and the interaction of the plasma with the wall. As a
consequence, less computationally demanding tools such as fluid simulations
(see, e.g., Refs. Stegmeir  _et al._ , 2018; De Oliveira  _et al._ , 2022;
Giacomin  _et al._ , 2022) based on the drift-reduced Braginskii model Zeiler,
Drake, and Rogers (1997), are often used to simulate the plasma dynamics in
the boundary. However, the Braginskii fluid approach remains limited to the
highly collisional region of the boundary, namely the scrape-off layer (SOL)
and cannot describe kinetic effects at lower collisionality.

To tackle the challenges of the boundary region, an approach is formulated in
Ref. Frei, Jorge, and Ricci, 2020 based on the Hermite-Laguerre expansion of
the full (full-F) distribution function, which leads to the evolution of an
arbitrary number of expansion coefficients, referred to as the gyro-moments
GMs. This approach features kinetic effects Jorge  _et al._ (2019); Frei,
Hoffmann, and Ricci (2022); Frei  _et al._ (2023), which are absent in
Braginskii-like fluid models, and collisional effects modeled using advanced
collision operators Jorge, Frei, and Ricci (2019); Frei  _et al._ (2021);
Frei, Ernst, and Ricci (2022); Frei  _et al._ (2023); Hoffmann, Frei, and
Ricci (2023a) with an arbitrary number of GMs. The ability to adjust the
number of GMs with a simple closure by truncation removes the need for ad-hoc
kinetic closures (Hammett and Perkins, 1990; Beer and Hammett, 1996) since the
accuracy of the model can be improved by increasing arbitrarily the number of
GMs and, therefore, the velocity-space resolution. So far, investigations
based on the GM approach are limited to the $\delta f$ regime (where only the
a priori small deviation of the distribution function from thermal equilibrium
is evolved (Brizard and Hahm, 2007)) showing convergence with a low number of
GMs, in particular at high collisionality Frei  _et al._ (2023). We notice
that the first nonlinear $\delta f$ simulations using advanced linearized
collision operators have been recently performed Hoffmann, Frei, and Ricci
(2023a, b) and that a similar approach has also been implemented to perform
$\delta f$ turbulent calculations focusing on the core region Mandell,
Dorland, and Landreman (2018); Mandell  _et al._ (2022). In this work, we
present the first full-F turbulent results that is based on the GM approach
using a flexible number of moments. In particular, we focus on simulations of
plasma turbulence in a linear plasma device.

Linear plasma devices, such as LAPD Gekelman  _et al._ (1991), HelCat Gilmore
_et al._ (2015) and RAID Furno  _et al._ (2017), are experiments that allow
for the investigation of basic plasma phenomena in a simplified magnetic
geometry characterized by the absence of magnetic gradients, curvature, and
shear Gekelman  _et al._ (1991); Tynan  _et al._ (2004); Castellanos  _et al._
(2005); Rogers and Ricci (2010); Furno  _et al._ (2017). Despite their
simplicity and the lack of kinetic effects such as trapped electrons, linear
plasma devices share some of the most important physical processes that occur
in the boundary of magnetic confinement devices. In fact, similar to the
boundary, the turbulent dynamics in a linear plasma device result from the
interplay of cross-field transport, parallel flows to the magnetic field, and
plasma losses at the end plates where a sheath forms due to plasma-wall
interactions. At the same time, the straight magnetic field lines in these
devices facilitate the development of new modelling tools, compared to complex
magnetic geometry characterizing the boundary of fusion devices. The modelling
in these devices is also simplified by the perpendicular incidence of the
magnetic field lines to the wall of the machine, which simplifies the sheath
boundary model compared to an oblique incidence Chodura (1982); Loizu  _et
al._ (2012); Geraldini, Parra, and Militello (2017) and by the low plasma
temperatures comparable to typical SOL values (e.g., $T_{i}\lesssim T_{e}\sim
6$ eV in typical LAPD discharges Friedman  _et al._ (2012)), which are ideal
for applying the full-F GM approach.

The low plasma temperature allows for a direct comparison of the GM approach
with validated fluid simulations Rogers and Ricci (2010); Popovich  _et al._
(2010); Fisher  _et al._ (2015); Ross  _et al._ (2019), which have been
considered because of the collisional conditions often met in, e.g., LAPD
experiments. Moreover, we note that, in addition to fluid simulations, LAPD
configuration was also chosen to perform the first full-F GK simulation in
open field lines with the `Gkeyll` code that uses a discontinuous-Galerkin
approach to discretize the velocity-space Shi  _et al._ (2017). LAPD turbulent
simulations using the GK `GENE` code are also reported in Ref. Pan  _et al._ ,
2018 based on the same physical model. In both cases, the GK model used for
the LAPD simulations neglects finite Larmor radius (FLR) effects and can,
therefore, be considered equivalent to a drift-kinetic model. As a matter of
fact, previous fluid and GK simulations of LAPD provide a comparison that make
linear plasma devices an ideal testbed to perform the first full-F turbulent
simulations using the GM approach with an arbitrary number of GMs.

In this work, we consider a simplified version of the full-F GM model derived
in Ref. Frei, Jorge, and Ricci, 2020. Consistently with Refs. Shi  _et al._ ,
2017 and Pan  _et al._ , 2018, we focus on the long-wavelength electrostatic
limit of the GK model to describe the ion dynamics, with ion-ion collisions
modeled using a simple nonlinear Dougherty Dougherty (1964) collision operator
(similar to the one used in Refs. Shi  _et al._ , 2017; Pan  _et al._ , 2018).
In addition, we neglect FLR effects and we refer to our ion model as a GK
model for consistency with previous works Shi  _et al._ (2017); Pan  _et al._
(2018), but emphasize its similarity to a drift-kinetic ion model. On the
other hand, electrons are assumed collisional, such that their dynamics can be
approximated by the drift-reduced Braginskii model Zeiler, Drake, and Rogers
(1997). This is in contrast with respect to previous GK simulations where
electrons are treated using the long-wavelength (without FLR effects) GK
equation. Hence, our model can be considered as a hybrid kinetic-fluid model,
with a GK ion and fluid electron description.

In contrast to previous GK simulations of linear devices Shi  _et al._ (2017);
Pan  _et al._ (2018), the ion GK equation is solved within the GM approach
where the full ion distribution function $F_{i}$ is expanded on a Hermite and
Laguerre polynomial basis. A parallel (to the magnetic field) velocity-space
coordinate shifted by the local ion parallel fluid velocity and the adiabatic
invariant are used to describe efficiently sonic ion parallel flows near the
end plates where the sheath forms. A full-F ion GM hierarchy equation for the
expansion coefficients is then derived. The ion full-F GM hierarchy equation
and the fluid electron model are coupled through a vorticity equation where
the Boussinesq approximation is considered. To incorporate the losses at the
end plates, Bohm sheath boundary conditions Loizu  _et al._ (2012) are
implemented in the parallel direction, which are equivalent to the ones used
in the previous Braginskii simulation of LAPD Rogers and Ricci (2010); Fisher
_et al._ (2015). Nonlinear simulations of LAPD are then performed with various
numbers of GMs. For comparison, a set of nonlinear turbulent simulations are
also performed using the two-fluid drift-reduced Braginskii equations Zeiler,
Drake, and Rogers (1997) (or simply Braginskii model), similarly to Refs.
Rogers and Ricci, 2010; Fisher  _et al._ , 2015, and using a reduced cold-ion
model derived from the full-F ion GM hierarchy.

The present results demonstrate that the full-F GM approach properly describes
fluctuations in an open-field line geometry. In fact, a detailed analysis
shows that turbulence, driven by a long perpendicular wavelength Kelvin-
Helmholtz (KH) instability, is in qualitative agreement with the Braginskii
model. The importance of the KH instability in determining the radial
turbulent transport in LAPD, as pointed out by Ref. Rogers and Ricci, 2010,
motivate the long-wavelength approximation and neglecting FLR effects in the
present work. Our findings exhibit weak dependence on the number of GMs
employed in the simulations and on the collisional regime. This is primarily
due to the absence of significant kinetic effects in LAPD and the fluid nature
of the KH instability, which governs radial turbulent transport and manifests
itself at long perpendicular wavelengths. The analysis of the velocity-space
representation of the ion distribution function demonstrates that the
amplitude of the GMs decays rapidly with the order of the polynomial when
collisions are considered. The present investigation also reveals that a
simple closure based on the truncation of the GM hierarchy is sufficient in
our case and has little effect on turbulence. It is important to note that the
purpose of these simulations is not to achieve a highly-fidelity and realistic
description of LAPD turbulence, but rather to establish confidence in the
applicability of the GM approach in full-F turbulent calculations with a
flexible number of GMs. Furthermore, a direct comparison with LAPD
experimental data Carter (2006); Fisher  _et al._ (2015) and with previous GK
simulations Shi  _et al._ (2017); Pan  _et al._ (2018) falls outside the scope
of our study, but will be addressed in future work.

We remark that the GM approach of Ref. Frei, Jorge, and Ricci, 2020 shares
some similarities with previous (full-F or $\delta f$) gyrofluid models
Brizard (1992); Beer and Hammett (1996); Snyder and Hammett (2001); Strintzi
and Scott (2004); Madsen (2013); Held  _et al._ (2016); Wiesenberger  _et al._
(2019); Galassi  _et al._ (2022); Wiesenberger and Held (2022, 2023), where a
fixed number of moments, usually up to the third-order velocity-space moments,
are evolved as dynamical variables. These models can be considered as the
limit of the GM approach where a finite number of moments is evolved and
properly designed closures to mimic kinetic effects, such as Landau damping
and FLR effects. Hence, the GM approach used in this work is similar to a
full-F gyrofluid model without FLR effects when the number of moments is kept
constant. Finally, although our model considers a closure by truncation, FLR
closures at arbitrary perpendicular wavelengths Held, Wiesenberger, and Kendl
(2020); Held and Wiesenberger (2022); Wiesenberger and Held (2023) and
collisional closures Jorge, Ricci, and Loureiro (2017); Wiesenberger and Held
(2022) for full-F gyrofluid models can also be considered with an arbitrary
number of GMs.

The paper is structured as follows. In Sec. II, we derive the ion full-F GM
hierarchy equation in a straight magnetic field and introduce the electron
fluid model, as well as the two-fluid drift-reduced Braginskii model. The
numerical implementation of the full-F GM hierarchy equation is detailed in
Sec. III. The results of the first full-F GM turbulent simulations are
presented in Sec. IV, which includes a detailed comparison with the Braginskii
simulations and an analysis of the ion distribution function. We conclude in
Sec. V. Appendix A reports on the derivation of the ion temperature equation
at high collisionality using the GM approach.

##  II Linear Plasma Device Model

In this section, we derive the full-F GM hierarchy equation for the ion
dynamics by expanding the ion distribution function onto a Hermite and
Laguerre polynomials basis. The hierarchy includes particle and energy sources
and a simple nonlinear long-wavelength Dougherty collision operator. The ion
full-F GM hierarchy is derived by neglecting FLR effects from gyro-average and
the difference between particle and gyrocenter fluid quantities, but includes
long-wavelength approximation of the polarization density in the
quasineutrality equation where the Boussinesq approximation is considered.
Ultimately, this allows us to remove the complexities of FLR effects in full-F
calculations using a flexible number of GMs. A reduced cold-ion model is also
considered, for comparison purposes, which is obtained analytically from the
full-F GM hierarchy in the cold-ion ($T_{i}=0$) limit. For the electron
dynamics, the Braginskii fluid equations are used to evolve the electron
density $n_{e}$, parallel velocity $U_{\parallel e}$, and temperature $T_{e}$.
A vorticity equation is derived using the Boussinesq approximation for the
electrostatic potential $\phi$, which couples the ion and electron models.
Finally, we present the two-fluid Braginskii model.

The simple magnetic geometry in a linear plasma device allows us to introduce
a simple coordinate system. In particular, assuming a rectangular shape of the
linear device cross-section, we define the cartesian coordinate system
$(x,y,z)$, such that the $(x,y)$ coordinates describe the plane perpendicular
to $\bm{B}$, while $z$ is the coordinate along the magnetic field lines. The
height and width of the perpendicular cross-section are $L_{x}$ and $L_{y}$,
respectively, and the length of the linear plasma device is $L_{z}$. The
constant in time and uniform in space magnetic field can simply be written as
$\bm{B}=\gradient\times\bm{A}=B\bm{z}$, where $\bm{A}$ is the magnetic vector
potential, which is constant in time, and $\bm{z}$ is the unit vector pointing
along the axis of the linear device.

This section is structured as follows. We describe the ion full-F model in
Sec. II.1 and we derive the ion full-F GM hierarchy equation in Sec. II.2. A
presentation of the reduced cold-ion model is then obtained from the GM
hierarchy equation in Sec. II.3. The fluid electron model follows in Sec. II.4
and the vorticity equation is derived in Sec. II.5. Sec. II.6 describes the
two-fluid Braginskii model used for comparison purposes and, finally, Sec.
II.7 details the Bohm sheath boundary conditions we use in our simulations.

###  II.1 Ion full-F model

Focusing on the electrostatic and neglecting FLR effects with constant and
straight magnetic field lines, the ion one-form $\Gamma_{i}$ Frei, Jorge, and
Ricci (2020) expressed in the gyrocenter coordinates
$\bm{Z}=(\bm{R}_{g},\mu,v_{\parallel},\theta)$, where $\bm{R}_{g}$ is the
gyrocenter position, $\mu=m_{i}v_{\perp}^{2}/(2B)$ is the magnetic moment and
$v_{\parallel}=\bm{b}\cdot\bm{v}$ is the velocity parallel to the magnetic
field with $\bm{b}=\bm{B}/B$, reduces to

$\displaystyle\Gamma_{i}(\bm{R}_{g},\mu,v_{\parallel},t)=q_{i}\bm{A}_{i}^{*}\cdot\dot{\bm{R}}_{g}-\frac{\mu
B}{\Omega_{i}}\dot{\theta}-m_{i}v_{\parallel}^{2}/2-q_{i}\Phi_{i},$ (1)

with $q_{i}\Phi_{i}=q_{i}\phi+\mu B$ and
$q_{i}\bm{A}_{i}^{*}=q_{i}\bm{A}+m_{i}v_{\parallel}\bm{b}$. We remark that, in
Eq. (1), the electrostatic potential $\phi$ is evaluated at the gyrocenter
position, i.e. $\phi=\phi(\bm{R}_{g})$, such that ion FLR effects are
neglected. From Eq. (1), we deduce the ion equations of motion,

$\displaystyle\dot{\bm{R}_{g}}$
$\displaystyle=\bm{b}v_{\parallel}+\frac{\bm{E}\times\bm{b}}{B},$ (2a)
$\displaystyle\dot{v}_{\parallel}$
$\displaystyle=\frac{q_{i}}{m_{i}}\bm{b}\cdot\bm{E},$ (2b)

and $\dot{\mu}=0$, with $v_{\parallel}=\dot{\bm{R}_{g}}\cdot\bm{b}$ and
$\bm{E}=-\gradient\phi$ the electric field. Eq. (2a) describes the parallel
streaming along the magnetic field lines and the perpendicular drift due to
the $\bm{E}\times\bm{B}$ velocity, while Eq. (2b) represents the acceleration
in the parallel direction associated with the electric field $\bm{E}$. Using
the equations of motion given in Eq. (2), the evolution equation of the full-F
(gyrophase-independent) ion distribution function,
$F_{i}=F_{i}(\bm{R}_{g},\mu,v_{\parallel},t)$, in the electrostatic limit of
the GK ion Boltzmann equation Brizard and Hahm (2007) is given by

$\displaystyle\frac{\partial}{\partial
t}\left(\mathcal{J}_{i}F_{i}\right)+\gradient\cdot\left(\mathcal{J}_{i}\dot{\bm{R}}_{g}F_{i}\right)+\frac{\partial}{\partial
v_{\parallel}}\left(\mathcal{J}_{i}F_{i}\dot{v}_{\parallel}\right)=\mathcal{J}_{i}\mathcal{C}_{i}+\mathcal{J}_{i}S_{i},$
(3)

where $\mathcal{J}_{i}=B/m_{i}$ is the gyrocenter phase-space Jacobian, which
is a constant in the case of linear devices. On the right-hand side of Eq.
(3),
$S_{i}=S_{i}(\bm{R}_{g},\mu,v_{\parallel})=S_{\mathcal{N}}+S_{\mathcal{E}}$
models the particle ($S_{N}$) and energy ($S_{E}$) sources and are defined by
Goerler  _et al._ (2011); Sarazin  _et al._ (2010)

$\displaystyle S_{\mathcal{N}}$
$\displaystyle=\mathcal{A}_{\mathcal{N}}F_{Mi},$ (4a) $\displaystyle
S_{\mathcal{E}}$ $\displaystyle=\mathcal{A}_{\mathcal{E}}\left(s_{\parallel
i}^{2}+x_{i}-\frac{3}{2}\right)F_{Mi},$ (4b)

respectively. We note that the sources, given in Eq. (4), neglect FLR and
polarization effects associated with the transformation from particle to
gyrocenter coordinate Wiesenberger and Held (2022). In Eq. (4), the
perpendicular normalized velocity-space coordinate is defined by $x_{i}=\mu
B/T_{i0}$, while we use $s_{\parallel i}=(v_{\parallel}-U_{\parallel
i})/v_{Ti}$ (with $v_{Ti}^{2}=2T_{i0}/m_{i}$ and $U_{\parallel
i}=\bm{b}\cdot\bm{u}_{i}=\int d\bm{v}F_{i}v_{\parallel}/N_{i}$ the ion
parallel fluid velocity) for the normalized parallel velocity-space
coordinate. Here, $T_{i0}$ is the reference ion temperature, which is assumed
constant in time and space. We remark that in previous fluid investigations of
LAPD, the low ion temperature assumption ($T_{i}\ll T_{e}$) is used and the
ion energy source is neglected.

The functions $\mathcal{A}_{\mathcal{N}}=\mathcal{A}_{\mathcal{N}}(x,y)$ and
$\mathcal{A}_{\mathcal{E}}=\mathcal{A}_{\mathcal{E}}(x,y)$ in Eq. (4) are
chosen to describe the spatial localization of the sources resulting from
ionization processes due to fast electrons and ions Tripathi, Pribyl, and
Gekelman (2011). Neglecting the presence of localized sources (not uniform in
$z$) near the end plates, we assume that these sources are a uniform in $z$
and radially localized with a top-hat-like shape. For instance,
$\mathcal{A}_{\mathcal{N}}(x,y)$ is given by Rogers and Ricci (2010)

$\displaystyle\mathcal{A}_{\mathcal{N}}(x,y)=\mathcal{A}_{\mathcal{N}0}0.5\left[1-\tanh\left(\frac{r-r_{s}}{L_{s}}\right)\right]+\mathcal{A}_{\mathcal{N}\infty},$
(5)

where $r=\sqrt{x^{2}+y^{2}}$ is the perpendicular distance from the center of
the device ($r=0$), $r_{s}$ is the radial extent of the plasma source,
$L_{s}>0$ is its typical source decay scale length, $\mathcal{A}_{N0}$ is a
positive and constant coefficient, which represents the physical particle
fuelling rate near the center of the device, while
$\mathcal{A}_{\mathcal{N}\infty}$ represents a small
($\mathcal{A}_{\mathcal{N}\infty}\ll\mathcal{A}_{\mathcal{N}0}$) positive and
constant particle source away from $r\sim r_{s}$ added for numerical reasons,
in particular, to avoid regions of small or negative plasma density. Similar
definitions for $\mathcal{A}_{\mathcal{E}}(x,y)$, $\mathcal{A}_{\mathcal{E}0}$
and $\mathcal{A}_{\mathcal{E}\infty}$ are used. In Eq. (4), we also introduce
a shifted Maxwellian distribution function defined by Brizard (1992); Madsen
(2013)

$\displaystyle F_{Mi}=\frac{N_{i}}{\pi^{3/2}v_{Ti}^{3}}e^{-s_{\parallel
i}^{2}}e^{-x_{i}},$ (6)

which assumes isotropic ion temperatures. We remark that the sources given in
Eq. (4) neglect non-Maxwellian sources, which can be relevant to fast ion
injections, and that the functional form of Eq. (4) has a straightforward
projection in the Hermite-Laguerre basis.

Finally, the term $\mathcal{C}_{i}$ in Eq. (3) is a full-F and nonlinear
collision operator model describing ion-ion collisions. In particular, we use
a long-wavelength (neglecting FLR effects) Dougherty collision operator
Dougherty (1964), given by

$\mathcal{C}_{i}=\nu_{i}\frac{\partial}{\partial\bm{v}}\cdot\left[\frac{2T_{i}}{m_{i}}\frac{\partial}{\partial\bm{v}}{F_{i}}-(\bm{v}-\bm{u}_{i})F_{i}\right],$
(7)

where $T_{i}=\int d\bm{v}F_{i}m_{i}(\bm{v}-\bm{u}_{i})^{2}/(3N_{i})$ and
$\bm{u}_{i}=\int d\bm{v}F_{i}\bm{v}/N_{i}$ are the ion temperature and mean
fluid velocity, and
$\nu_{i}=4\sqrt{\pi}N_{i}q_{i}^{4}\ln\lambda/(3m_{i}^{1/2}T_{i0}^{3/2})$ is
the ion-ion collision frequency, which is constant in the present work. The
effects of ion-electron collisions are neglected in Eq. (3) since they occur
on a time scale larger by, at least, a factor proportional to
$\sqrt{m_{i}/m_{e}}$ than ion-ion collisions.

Eq. (3) is equivalent to the ion GK model used in previous GK turbulent
simulations of LAPD, implemented in the `Gkeyll` Shi  _et al._ (2017) and in
the `GENE` Pan  _et al._ (2018) codes. Both implementations use the same
nonlinear Dougherty collision operator for ion-ion collisions (given in Eq.
(7)) and neglect ion-electron collisions. The `Gkeyll` code employs a
discontinuous-Galerkin approach and `GENE` uses a finite-volume method to
discretize the velocity-space coordinates $(v_{\parallel},\mu)$, while our
work uses the GM approach to simulate the full-F ion distribution function
$F_{i}$. To our knowledge, this is the first time such a moment approach is
applied to perform nonlinear full-F turbulent simulations with an arbitrary
number of ion GMs.

###  II.2 Full-F ion GM Hierarchy Equation

Following Ref. Frei, Jorge, and Ricci, 2020, we perform the GM expansion of
the full-F ion distribution function, $F_{i}$. More precisely, we expand
$F_{i}$ onto a set of Hermite, $H_{p}(s_{\parallel i})$, and Laguerre,
$L_{j}(x_{i})$, velocity-space polynomials Gradshteyn and Ryzhik (2014), such
that

$\displaystyle
F_{i}=\sum_{p=0}^{\infty}\sum_{j=0}^{\infty}\mathcal{N}^{pj}\frac{H_{p}(s_{\parallel
i})L_{j}(x_{i})}{\sqrt{2^{p}p!}}\frac{F_{Mi}}{N_{i}},$ (8)

where $\mathcal{N}^{pj}$ are the ion GMs, evaluated by using the Hermite and
Laguerre orthogonality relations Gradshteyn and Ryzhik (2014)

$\displaystyle\mathcal{N}^{pj}=2\pi\int_{-\infty}^{\infty}dv_{\parallel}\int_{0}^{\infty}d\mu\frac{B}{m_{i}}F_{i}\frac{H_{p}(s_{\parallel
i})L_{j}(x_{i})}{\sqrt{2^{p}p!}}.$ (9)

By introducing the GM projector Madsen (2013); Jorge, Ricci, and Loureiro
(2017); Mandell, Dorland, and Landreman (2018); Frei, Jorge, and Ricci (2020)

$\displaystyle\norm{\chi}^{pj}=2\pi\int_{-\infty}^{\infty}dv_{\parallel}\int_{0}^{\infty}d\mu\chi\frac{B}{m_{i}}F_{i}\frac{H_{p}(s_{\parallel
i})L_{j}(x_{i})}{\sqrt{2^{p}p!}},$ (10)

with $\chi=\chi(\bm{R}_{g},\mu,v_{\parallel},t)$ being an arbitrary gyrocenter
phase-space function, we find $\mathcal{N}^{pj}=\norm{1}^{pj}$ from Eq. (10).
We remark that, in Eq. (8), the shift of the parallel velocity coordinate
$s_{\parallel i}$, appearing in $F_{Mi}$ defined in Eq. (6) and in the
argument of the Hermite polynomial $H_{p}$, is necessary to ensure good
convergence property of the GM approach with respect to the number of GMs in
Eq. (8), in particular in the presence of sonic ion flows (see Sec. IV.3).
These flows appear at the sheath entrance where ions are accelerated to the
ion sound speed (see Sec. II.7). Additionally, we note that $F_{Mi}$, defined
in Eq. (6), is assumed to have the same parallel and perpendicular
temperature, $T_{\parallel i}=T_{\perp i}=T_{i0}$. The assumption of an
isotropic Maxwellian distribution function in Eq. (8) is justified by the
large ion-ion collision frequency typically found in a linear plasma device
(where $T_{i}\lesssim 1$ eV) compared to the boundary region in fusion devices
(where $T_{i}\gtrsim 10$ eV). The absence of strong external energy sources
driving temperature anisotropy in LAPD experiments supports this assumption
(see Eq. (4)).

The lowest-order GMs can be related to fluid ion gyrocenter quantities, such
as the ion gyrocenter density $N_{i}$, the ion parallel velocity $U_{\parallel
i}$, and the ion parallel and perpendicular pressure and temperature
$P_{\parallel i}=T_{\parallel i}N_{i}$ and $P_{\perp i}=T_{\perp i}N_{i}$,
respectively. Indeed, using Eq. (9), we derive that

$\displaystyle N_{i}$ $\displaystyle=\mathcal{N}^{00},$ (11a)
$\displaystyle\mathcal{N}^{10}$ $\displaystyle=0,$ (11b) $\displaystyle
P_{\parallel i}=N_{i}T_{\parallel i}$
$\displaystyle=2\pi\int_{-\infty}^{\infty}dv_{\parallel}\int_{0}^{\infty}d\mu
BF_{i}(v_{\parallel}-U_{\parallel i})^{2}$
$\displaystyle=T_{i0}\left(\sqrt{2}\mathcal{N}^{20}+N_{i}\right),$ (11c)
$\displaystyle P_{\perp i}=N_{i}T_{\perp i}$
$\displaystyle=2\pi\int_{-\infty}^{\infty}dv_{\parallel}\int_{0}^{\infty}d\mu\frac{B}{m_{i}}\mu
BF_{i}$ $\displaystyle=T_{i0}(N_{i}-\mathcal{N}^{01}),$ (11d)

with the total ion temperature defined by $T_{i}=\int
d\bm{v}F_{i}m_{i}(\bm{v}-\bm{b}U_{\parallel i})^{2}/(3N_{i})=(T_{\parallel
i}+2T_{\perp
i})/3=T_{i0}(\sqrt{2}\mathcal{N}^{20}+3N_{i}-2\mathcal{N}^{01})/(3N_{i})$. We
remark that Eq. (11b) is a direct consequence of our choice of using a shifted
parallel velocity-space coordinate $s_{\parallel i}$ in Eq. (8). In Appendix
A, we use the definitions of Eq. (11) to derive the Braginskii ion temperature
equation. We note that the ion gyrocenter fluid quantities given in Eq. (11)
are equivalent to the particle ones since FLR and polarization effects Jorge,
Ricci, and Loureiro (2017); Jorge, Frei, and Ricci (2019); Wiesenberger and
Held (2022), associated with the transformation from particle to gyrocenter
coordinates, are neglected in our model.

We now derive the full-F GM hierarchy equation describing the evolution of an
arbitrary number of GMs, $\mathcal{N}^{pj}$. This is obtained by projecting
the ion full-F equation given in Eq. (3) onto the Hermite-Laguerre basis. In
addition, we normalize time $t$ to $R/c_{s0}$ (with
$c_{s0}=\sqrt{T_{e0}/m_{i}}$ the ion sound speed evaluated at the reference
constant electron temperature $T_{e0}$ and $R$ the radial extension of the
plasma chamber in the direction perpendicular to $\bm{B}$), the potential
$\phi$ to $T_{e0}/e$, the parallel and perpendicular spatial scales to $R$ and
$\rho_{s0}=c_{s0}/\Omega_{i}$, respectively. We also normalize the ion and
electron densities, $N_{i}$ and $N_{e}$, to the constant reference density
$N_{0}$, the parallel electron velocity $U_{\parallel e}$ to $c_{s0}$, and the
electron temperature, $T_{e}$, to $T_{e0}$. In addition, we assume $q_{i}=+e$,
considering a hydrogen plasma. Hence, we derive the normalized ion GM
hierarchy equation, which describes the evolution of the GMs
$\mathcal{N}^{pj}$, i.e.

$\displaystyle\frac{\partial}{\partial
t}\mathcal{N}^{pj}+\sqrt{\frac{p}{\tau_{i}}}\mathcal{N}^{p-1j}\frac{\partial}{\partial
t}U_{\parallel i}+\gradient\cdot\norm{\dot{\bm{R}}_{g}}^{pj}$
$\displaystyle+\sqrt{\frac{p}{\tau_{i}}}\norm{\dot{\bm{R}_{g}}}^{p-1j}\cdot\gradient
U_{\parallel
i}-\sqrt{\frac{p}{\tau_{i}}}\norm{\dot{v}_{\parallel}}^{p-1j}=\mathcal{C}^{pj}_{i}+S_{\mathcal{N}}^{pj}+S_{\mathcal{E}}^{pj},$
(12)

where the GM projections are given by

$\displaystyle\gradient\cdot\norm{\dot{\bm{R}_{g}}}^{pj}$
$\displaystyle=\sqrt{\tau_{i}}\partial_{z}\left(\sqrt{p+1}\mathcal{N}^{p+1j}+\sqrt{p}\mathcal{N}^{p-1j}\right)$
$\displaystyle+\partial_{z}\left(U_{\parallel
i}\mathcal{N}^{pj}\right)+\frac{1}{\rho_{*}}\left[\phi,\mathcal{N}^{pj}\right],$
(13a)
$\displaystyle\sqrt{\frac{p}{\tau_{i}}}\norm{\dot{\bm{R}_{g}}}^{p-1j}\cdot\gradient
U_{\parallel
i}=\left(p\mathcal{N}^{pj}+\sqrt{p(p-1)}\mathcal{N}^{p-2j}\right.$
$\displaystyle\left.+\sqrt{\frac{p}{\tau_{i}}}U_{\parallel
i}\mathcal{N}^{p-1j}\right)\partial_{z}U_{\parallel
i}+\sqrt{\frac{p}{\tau_{i}}}\frac{1}{\rho_{*}}\mathcal{N}^{p-1j}\left[\phi,U_{\parallel
i}\right],$ (13b)
$\displaystyle\sqrt{\frac{p}{\tau_{i}}}\norm{\dot{v}_{\parallel}}^{p-1j}$
$\displaystyle=-\mathcal{N}^{p-1j}\sqrt{\frac{p}{\tau_{i}}}\partial_{z}\phi,$
(13c)

with $\rho^{*}=\rho_{s0}/R$ and $\tau_{i}=T_{i0}/T_{e0}$. In Eq. (13), we
introduce the Poisson bracket operator that is
$\left[f,g\right]=\partial_{x}f\partial_{y}g-\partial_{y}f\partial_{x}g$. The
GM expansions of the particle and energy sources, $S_{N}^{pj}$ and
$S_{\mathcal{E}}^{pj}$, are given by

$\displaystyle
S_{\mathcal{N}}^{pj}=\mathcal{A}_{\mathcal{N}}\delta_{p}^{0}\delta_{j}^{0}=S_{\mathcal{N}},$
(14a) $\displaystyle
S_{\mathcal{E}}^{pj}=\mathcal{A}_{\mathcal{E}}\left(\frac{\delta_{p}^{2}\delta_{j}^{0}}{\sqrt{2}}-\delta_{p}^{0}\delta_{j}^{1}\right),$
(14b)

respectively. The contribution of Eq. (14a) to the ion density equation
obtained with $(p,j)=(0,0)$ represents particle sources, while the non-
vanishing terms with $(p,j)=(2,0)$ and $(0,1)$ in Eq. (14b) are associated
with parallel and perpendicular energy sources.

Finally, we express the nonlinear Dougherty collision operator in terms of
GMs. We first express Eq. (7) in terms of the velocity-space coordinates
$(s_{\parallel i},x_{i})$ and project it onto the Hermite-Laguerre basis. This
yields

$\displaystyle\mathcal{C}^{pj}_{i}$
$\displaystyle=\nu_{i}\left[-(p+2j)\mathcal{N}^{pj}+\left(T_{i}-1\right)\right.$
$\displaystyle\left.\times\left(\sqrt{p(p-1)}\mathcal{N}^{p-2j}-2j\mathcal{N}^{pj-1}\right)\right],$
(15)

where $T_{i}$ is expressed in terms of the GMs using Eq. (11). The nonlinear
Dougherty collision operator conserves particles ($\mathcal{C}_{i}^{00}=0$),
momentum ($\mathcal{C}_{i}^{10}=0$) and energy
($\mathcal{C}_{i}^{20}=\sqrt{2}\mathcal{C}_{i}^{01}$). While simpler in form
compared to the GM expansion of the nonlinear Fokker-Planck Landau collision
operator Jorge, Frei, and Ricci (2019), the Dougherty collision operator
constitutes an initial step to incorporate advanced collisional effects in the
nonlinear and full-F ion GM hierarchy equation. The numerical implementation
of the nonlinear Fokker-Planck Landau collision operator Jorge, Frei, and
Ricci (2019) will be considered in future work.

To obtain the time evolution of the GMs $\mathcal{N}^{pj}$, it is necessary to
derive an explicit expression for the time derivative of the ion parallel
velocity, $\partial_{t}U_{\parallel i}$ which appears in Eq. (II.2) and
resulting from the use of the shifted parallel velocity-space coordinate
$s_{\parallel i}$. By setting $(p,j)=(1,0)$ in Eq. (II.2) and using the fact
that $\mathcal{N}^{10}$ vanishes exactly (see Eq. (11b)), we derive the
desired expression for $\partial_{t}U_{\parallel i}$ given by

$\displaystyle N_{i}\partial_{t}U_{\parallel
i}+\frac{N_{i}}{\rho_{*}}\left[\phi,U_{\parallel
i}\right]+\tau_{i}\partial_{z}P_{\parallel i}+N_{i}U_{\parallel
i}\partial_{z}U_{\parallel i}+N_{i}\partial_{z}\phi=0,$ (16)

where the parallel ion pressure $P_{\parallel i}$ is expressed in terms of GMs
according to Eq. (11c).

We note that the full-F GM hierarchy equation, given in Eq. (II.2), can also
be derived from the electromagnetic full-F GM hierarchy equation described in
Ref. Frei, Jorge, and Ricci, 2020. This is achieved by considering the
electrostatic limit, neglecting FLR effects, and assuming isotropic and
constant ion temperatures in $F_{Mi}$. We remark that, if $F_{Mi}$ is
normalized with the local parallel and perpendicular ion temperatures
($T_{\parallel i}$ and $T_{\perp i}$ in Eq. (11c) and Eq. (11d), respectively)
as in Ref. Frei, Jorge, and Ricci, 2020, additional terms proportional to
their time derivatives and gradients appear in Eq. (II.2). These terms might
be important in the case of large temperature variations such as the ones in
the boundary, but can be neglected in the case of LAPD.

Notably, the GMs with different $p$ are coupled in Eq. (II.2) due to the
parallel streaming terms, associated with the ion Landau damping. On the other
hand, the GMs with different $j$ are only coupled through the collision
operator (see Eq. (II.2)) as our model neglects FLR effects yielding
additional coupling in $j$ Frei, Hoffmann, and Ricci (2022). As a result, a
few Laguerre GMs are expected to be sufficient in our nonlinear turbulent
simulations.

To carry out the numerical turbulent simulations presented here, a simple
closure by truncation is applied to the GM hierarchy equation. More precisely,
we set $\mathcal{N}^{pj}=0$ for all $(p,j)>(P,J)$ with $0\leq P,J<\infty$. The
full-F GM hierarchy equation enables us to perform turbulent simulations of
LAPD using an arbitrary number of GMs. Different values of $(P,J)$ are
considered in Sec. IV where we demonstrate that the closure by truncation is
sufficient to perform full-F turbulent simulations in our case.

###  II.3 Cold-ion reduced model

We consider here the cold-ion limit of the full-F GM hierarchy and derive a
simplified model, similar to the one used in previous turbulent investigations
of linear devices based on fluid models (see, e.g., Refs. Rogers and Ricci,
2010; Popovich  _et al._ , 2010; Fisher  _et al._ , 2015) where the effects of
finite ion temperature $T_{i}$ are neglected.

In the cold ion limit ($T_{i}=0$), only the GM $\mathcal{N}_{i}^{00}$ and
$U_{\parallel i}$, associated with the ion gyrocenter density and the parallel
ion velocity, need to be evolved and the contribution from the parallel ion
pressure $P_{\parallel i}$ in Eq. (16) can be neglected. As a consequence, the
ion GM hierarchy equation given in Eq. (II.2) reduces to the ion gyrocenter
continuity equation for $N_{i}$ and to the ion parallel momentum equation for
$U_{\parallel i}$ Wiesenberger and Held (2020, 2022); Galassi  _et al._
(2022), i.e.

$\displaystyle\frac{\partial}{\partial
t}N_{i}+\frac{1}{\rho^{*}}\left[\phi,N_{i}\right]+\partial_{z}\left(U_{\parallel
i}N_{i}\right)=S_{\mathcal{N}},$ (17a) $\displaystyle\partial_{t}U_{\parallel
i}+\frac{1}{\rho_{*}}\left[\phi,U_{\parallel i}\right]+U_{\parallel
i}\partial_{z}U_{\parallel i}+\partial_{z}\phi=0,$ (17b)

respectively. We remark that the particle and momentum conservation of the
collision operator is used in deriving Eq. (17).

###  II.4 Electron fluid model

We use the Braginskii model to evolve the electron dynamics, avoiding the
evolution of their distribution function, in contrast to Refs. Shi  _et al._ ,
2017; Pan  _et al._ , 2018. The fluid approach for the electrons is justified
when the electron collision frequency is much larger than the ion collision
frequency and electron FLR effects are negligible for modes developing at
$k_{\perp}\rho_{s}\sim 1$, which is the case of LAPD experiments. In the
absence of electron FLR effects, the electron particle and gyrocenter position
coinicide such that no distinction is made between the electron particle and
gyrocenter fluid quantities. Hence, the time evolution of the electron density
$n_{e}$, electron parallel velocity $U_{\parallel e}$, and temperature $T_{e}$
is determined by the continuity equation, the generalized Ohm’s law, and the
temperature equation, respectively. These equations are given by

$\displaystyle\partial_{t}n_{e}$
$\displaystyle+\frac{1}{\rho^{*}}\left[\phi,n_{e}\right]+\partial_{z}\left(U_{\parallel
e}n_{e}\right)=S_{\mathcal{N}},$ (18a) $\displaystyle\partial_{t}U_{\parallel
e}$ $\displaystyle+\frac{1}{\rho^{*}}\left[\phi,U_{\parallel
e}\right]+U_{\parallel e}\partial_{z}U_{\parallel
e}=\frac{m_{i}}{m_{e}}\left[\nu_{\parallel}J_{\parallel}+\partial_{z}\phi\right.$
$\displaystyle\left.-\frac{T_{e}}{n_{e}}\partial_{z}n_{e}-1.71\partial_{z}T_{e}\right],$
(18b) $\displaystyle\partial_{t}T_{e}$
$\displaystyle+\frac{1}{\rho^{*}}\left[\phi,T_{e}\right]+U_{\parallel
e}\partial_{z}T_{e}=\frac{2}{3}T_{e}\left(\frac{0.71}{N_{e}}\partial_{z}J_{\parallel}-\partial_{z}U_{\parallel
e}\right)$ $\displaystyle+\partial_{z}\left(\chi_{\parallel
e}\partial_{z}T_{e}\right)+S_{T_{e}},$ (18c)

where the normalized parallel electrical resistivity and electron thermal
conductivity are given by $\nu_{\parallel}=\nu_{0}/T_{e}^{3/2}$ and
$\chi_{\parallel e}=1.075T_{e}^{5/2}/\nu_{0}$, respectively. Here,
$\nu_{0}=4\sqrt{2\pi}e^{4}n_{e0}R\sqrt{m_{e}}\ln\lambda/[3c_{s0}m_{i}T_{e0}^{3/2}1.96]$
is the normalized electron collisionality. On the right-hand side of Eq. (18a)
and Eq. (18c), $S_{N}$ and $S_{T_{e}}$ are the normalized density and
temperature sources. In Eq. (18), the parallel electrical current is
$J_{\parallel}=n_{e}(U_{\parallel i}-U_{\parallel e})$.

###  II.5 Vorticity equation

We now obtain the vorticity equation that governs the evolution of the
electrostatic potential $\phi$. This equation imposes the charge conservation
constraint to the time evolution of the plasma densities and electrical
currents.

To derive the vorticity equation, we consider the quasineutrality condition in
the long-wavelength limit, given by Brizard and Hahm (2007); Frei, Jorge, and
Ricci (2020)

$\displaystyle-
en_{e}+q_{i}N_{i}=-\gradient\cdot\left(\frac{q_{i}^{2}N_{i}}{m_{i}\Omega_{i}^{2}}\gradient_{\perp}\phi\right).$
(19)

We remark that Eq. (19) neglects FLR effects (associated with the gyro-
average) but retains the polarization effect in its right-hand side. Including
high-order FLR contributions to the gyro-average and polarization terms, which
are important for predicting of small-scale fluctuations Held and Wiesenberger
(2022), requires the development of FLR closures at arbitrary wavelength (see,
e.g., Ref. Held, Wiesenberger, and Kendl, 2020) using an arbitrary number of
GMs, a task left to future work. We also notice that Eq. (19) is equivalent to
the quasineutrality condition used in previous GK turbulent simulations of
LAPD Shi  _et al._ (2017); Pan  _et al._ (2018) if the Boussinesq
approximation is used, i.e. if $N_{i}$ is approximated by $N_{0}$ on the
right-hand side of Eq. (19). The Boussinesq approximation is widely used in
fluid codes Ross  _et al._ (2019); Giacomin  _et al._ (2022) and is also
considered below to derive the vorticity equation.

While Eq. (19) can be solved to obtain $\phi$ given the electron and ion
densities, $n_{e}$ and $N_{i}$ respectively, we use a vorticity equation
instead which is often considered in turbulent fluid codes Rogers and Ricci
(2010); Popovich  _et al._ (2010); Fisher  _et al._ (2015); Ross  _et al._
(2019). The vorticity equation is derived by taking the time derivative of the
quasineutrality equation given in Eq. (19) and by using the electron and ion
continuity equations, given in Eqs. (17a) and (18a), respectively.
Furthermore, we use the Boussinesq approximation, such that

$\displaystyle\gradient\cdot\left(\frac{q_{i}^{2}N_{i}}{m_{i}\Omega_{i}^{2}}\gradient_{\perp}\phi\right)\simeq\frac{q_{i}^{2}N_{i}}{m_{i}\Omega_{i}^{2}}\Omega,$
(20)

where we introduce the vorticity variable $\Omega=\gradient^{2}_{\perp}\phi$.
The vorticity equation is then

$\displaystyle-\partial_{t}\Omega-\frac{1}{\rho^{*}}\left[\phi,\Omega\right]-\partial_{z}(U_{\parallel
i}\Omega)+\frac{1}{N_{i}}\partial_{z}J_{\parallel}=0,$ (21)

We note that the effects of the Boussinesq approximation on plasma turbulence
are the subject of previous studies Yu, Krasheninnikov, and Guzdar (2006);
Bodi  _et al._ (2011); Paruta  _et al._ (2018); Ross  _et al._ (2019); Held
and Wiesenberger (2022). While it might not be justified in LAPD when steep
density gradients are present, it allows us to reduce the computational cost
of our simulations when inverting the two-dimensional Laplacian to obtain
$\phi$ from the vorticity variable $\Omega$. We use the vorticity equation
given in Eq. (21) to evolve $\Omega$ when considering the full-F ion GM
hierarchy equation and the cold ion models, given in Eqs. (II.2) and (17)
respectively, coupled to the fluid electron model in Eq. (18).

###  II.6 Two-fluid Braginskii fluid model

We finally introduce the two-fluid Braginskii fluid model Zeiler, Drake, and
Rogers (1997), valid in the high-collisional regime, for comparison with the
full-F ion GM hierarchy equation and the cold-ion model.

We first note that, in the two-fluid Braginskii model, quasineutrality is
assumed (neglecting the polarization term in Eq. (19)) and, as a consequence,
the electron particle density is used as an independent variable, such that
$n_{e}\simeq N_{i}$. On the other hand, ion polarization effects are retained
by the presence of the polarization drift (neglected in Eq. (2a)) in the ion
particle continuity equation, from which the vorticity equation (given below)
is derived for the electrostatic potential. Furthermore, since FLR effects
associated with the difference between particle and gyrocenter positions are
neglected in the two-fluid Braginskii fluid model, the ion gyrocenter fluid
quantities, such as ion parallel velocity $U_{\parallel i}$ and ion
temperature $T_{i}$, are assumed to be equal to particle fluid quantities.

In addition to the fluid electron fluid equations for $n_{e}$, $U_{\parallel
e}$ and $T_{e}$ already described in Sec. II.4, the two-fluid Braginskii
equations prescribe a parallel ion momentum equation to evolve $U_{\parallel
i}$, an ion temperature equation to evolve $T_{i}$, and vorticity equations
for $\Omega$. These equations are given by

$\displaystyle\partial_{t}U_{\parallel
i}+\frac{1}{\rho^{*}}\left[\phi,U_{\parallel i}\right]+U_{\parallel
i}\partial_{z}U_{\parallel i}$
$\displaystyle=-\partial_{z}T_{e}-\tau_{i}\partial_{z}T_{i}$
$\displaystyle-(T_{e}+\tau_{i}T_{i})\frac{\partial_{z}n_{e}}{n_{e}},$ (22a)
$\displaystyle\partial_{t}T_{i}+\frac{1}{\rho^{*}}\left[\phi,T_{i}\right]+U_{\parallel
i}\partial_{z}T_{i}$ $\displaystyle=+\frac{2}{3}T_{i}\left[\left(U_{\parallel
i}-U_{\parallel
e}\right)\frac{\partial_{z}n_{e}}{n_{e}}-\partial_{z}U_{\parallel e}\right]$
$\displaystyle+\partial_{z}\left(\chi_{\parallel
i}\partial_{z}T_{i}\right)+\frac{\mathcal{A}_{\mathcal{E}}}{n_{e}}+(1-T_{i})\frac{\mathcal{A}_{\mathcal{N}}}{n_{e}},$
(22b)
$\displaystyle\partial_{t}\Omega+\tau_{i}\partial_{t}\gradient_{\perp}^{2}T_{i}$
$\displaystyle=\frac{1}{n_{e}}\partial_{z}J_{\parallel}-\frac{1}{\rho^{*}}\left[\phi,\Omega+\tau_{i}\gradient_{\perp}^{2}T_{i}\right]$
$\displaystyle-U_{\parallel
i}\partial_{z}\left(\Omega+\tau_{i}\gradient_{\perp}^{2}T_{i}\right).$ (22c)

respectively. In Eq. (22b), $\chi_{\parallel
i}=1.32\sqrt{m_{e}/m_{i}}(\tau_{i}T_{i})^{5/2}/\nu_{0}$ is the normalized
parallel ion thermal conductivity. In Eq. (22b), the two last terms are the
ion temperature sources associated with the energy source $S_{\mathcal{E}}$
(see Eq. (14b)), which appears on the right-hand side of Eq. (3). In Appendix
A, we demonstrate analytically how these terms are determined by deriving the
ion temperature equation, Eq. (22b), from the ion GM hierarchy equation.

In contrast to the cold-ion model given in Eq. (17), the parallel electric
field $\partial_{z}\phi$, appearing in the parallel ion momentum equation, Eq.
(17b), is approximated in Eq. (22a) by the electron parallel pressure
gradient, such that $\partial_{z}\phi\simeq\partial_{z}P_{e}$ with
$P_{e}=n_{e}T_{e}$ (see Eq. (18b)). We also note that the terms proportional
to the Laplacian of the ion temperature, i.e.
$\tau_{i}\gradient_{\perp}^{2}T_{i}$, present in the vorticity equation, Eq.
(22), are absent in Eq. (21), which is derived from the GM hierarchy equation.
Indeed, these terms are associated with the long-wavelength FLR effect
correction, which are neglected in Eq. (21) as FLR effects are omitted. While
the ion temperature in LAPD experiments is generally lower than the electron
temperature, it has been shown that FLR effects can be important in this
regime. Further investigations are hence required to investigate the effects
of finite ion temperature and related FLR effects on small scale fluctuations
in LAPD.

###  II.7 Boundary conditions

Boundary conditions are required for the ion GMs, $\mathcal{N}^{pj}$, the
electron fluid quantities, $N_{e}$, $U_{\parallel e}$, $T_{e}$, and the
potential $\phi$ in the perpendicular $(x,y)$ plane at $x=\pm L_{x}/2$ and
$y=\pm L_{y}/2$ and at the end plates located in the $z$ direction at $z=\pm
L_{z}/2$, where a sheath forms due to the plasma-wall interaction.

At $x=\pm L_{x}/2$ and $y=\pm L_{y}/2$, homogenous Neumann boundary conditions
are used for all quantities. These ad-hoc boundary conditions have a
negligible effect on plasma turbulence near the center of the device as they
are imposed at a distance sufficiently large from the center of the device. On
the other hand, the boundary conditions in the $z$ direction have an important
impact since the formation of a Debye sheath is observed when the magnetic
field lines intercept the end plates that control the plasma losses Stangeby
(1995). Since the sheath region cannot be modeled by the field equations
derived in Sec. II.5 (the GK formalism is violated in this region), the sheath
is modeled in our simulations by a set of appropriate boundary conditions
imposed at the sheath entrance.

In previous GK simulations of LAPD Shi  _et al._ (2017); Pan  _et al._ (2018),
a conducting wall is considered. Accordingly, the fraction of electrons that
cross the sheath and are lost being absorbed by the walls is determined by the
value of the potential at the sheath entrance. This fraction is imposed by
evaluating the cutoff velocity of the electron distribution function
numerically. Leveraging the GM approach, we use the standard fluid Bohm
boundary conditions Stangeby (1995) which sets the value of the parallel
electron and ion velocities, $U_{\parallel e}$ and $U_{\parallel i}$, at the
sheath entrance. Therefore, we assume that Loizu  _et al._ (2012); Mosetto
_et al._ (2015)

$\displaystyle U_{\parallel e}(x,y,z=\pm L_{z}/2)$
$\displaystyle=\pm\sqrt{T_{e,s}}e^{\Lambda-\phi_{s}/T_{e,s}},$ (23a)
$\displaystyle U_{\parallel i}(x,y,z=\pm L_{z}/2)$ $\displaystyle=\pm c_{s}$
$\displaystyle=\pm\sqrt{T_{e,s}}\sqrt{1+\tau_{i}T_{i,s}/T_{e,s}},$ (23b)

with $\Lambda=\log m_{i}/(2m_{e})\simeq 3$ for hydrogen plasmas. In Eq. (23),
$T_{e,s}$ and $T_{i,s}$ are the electron and ion temperatures evaluated at the
sheath entrance, i.e. $T_{e,s}=T_{e}(x,y,z=\pm L_{z}/2)$ and
$T_{i,s}=T_{i}(x,y,z=\pm L_{z}/2)$, and, similarly $\phi_{s}=\phi(x,y,z=\pm
L_{z}/2)$. We notice that the boundary conditions in Eq. (23) reduce to the
ones used in Ref. Rogers and Ricci, 2010 when $T_{i}\ll T_{e}$ and correspond
to the ones used in SOL turbulent simulations using the drift-reduced
Braginskii model Giacomin  _et al._ (2022). For the remaining quantities, we
assume, for simplicity, that the gradients of electron density, $n_{e}$,
electron temperature, $T_{e}$, ion GMs, $\mathcal{N}^{pj}$, and electrostatic
potentials, $\phi$, vanish along the direction of the magnetic field at the
sheath entrance, i.e. homogenous Neumann boundary conditions are imposed at
$z=\pm L_{z}/2$. While the homogenous Neumann boundary conditions considered
here are sufficient to ensure the numerical stability of the present
simulations, further investigations are needed to develop first-principles
sheath boundary conditions for the GM approach. In particular, the analytical
procedure outlined in, e.g., Refs. Loizu  _et al._ , 2012; Mosetto  _et al._ ,
2015, can be extended to an arbitrary number of GMs and kinetic sheath
boundary conditions can also be developed (Riemann, 1981; Baalrud  _et al._ ,
2015; Churchill  _et al._ , 2017). Magnetic field lines intercept the machine
wall with a small oblique angle in fusion devices, further complicating the
treatment of the sheath boundary conditions Chodura (1982); Geraldini, Parra,
and Militello (2017).

##  III Numerical Implementation

To solve the full-F ion GM hierarchy given in Eq. (II.2) coupled with the
electron fluid model in Eq. (18) and the vorticity equation in Eq. (21), we
have developed a new three-dimensional full-F simulation code. This code
solves the turbulent dynamics for an arbitrary number of ion GMs. It also
solves the cold ion model (Sec. II.3) and the two-fluid Braginskii model (Sec.
II.6) for comparison with the GM results.

To evolve the plasma dynamics, we employ similar numerical algorithms as the
two-fluid `GBS` code Giacomin  _et al._ (2022). More precisely, an explicit
fourth-order Runge-Kutta time-stepping scheme is used. The perpendicular and
parallel directions are discretized using a uniform cartesian grid in the
$(x,y,z)$ coordinates with the $x$, $y$, and $z$ directions discretized using
$N_{x}$, $N_{y}$ and $N_{z}$ points uniformly distributed between the
intervals $[-L_{x}/2,+L_{x}/2]$, $[-L_{y}/2,+L_{y}/2]$ and
$[-L_{z}/2,L_{z}/2]$, respectively. The Poisson bracket operator,
$\left[f,g\right]=\bm{b}\times\gradient f\cdot\gradient
g=\partial_{x}f\partial_{y}g-\partial_{y}f\partial_{x}g$, with
$\bm{b}=\bm{B}/B=\bm{e}_{z}$, is evaluated by using a fourth-order Arakawa
method (Arakawa, 1997). The numerical evaluation of the other spatial
operators appearing in the GM hierarchy equation is based on a fourth-order
and centered finite difference scheme, resulting in a $5$-points centered
stencil Giacomin  _et al._ (2022). To avoid checkerboard patterns (Paruta  _et
al._ , 2018), the grid used to evolve the parallel velocities, $U_{\parallel
e}$ and $U_{\parallel i}$, and the GMs $\mathcal{N}^{pj}$ with odd $p$, is
staggered to the left along the $z$-direction by $\Delta z/2$ ($\Delta z$ is
the grid spacing) with respect to the grid where the other fluid quantities,
i.e. $n_{e}$, $T_{e}$, $\Omega$ (and thus $\phi$), and the GMs
$\mathcal{N}^{pj}$ with even $p$ are evaluated. Fourth-order interpolation
techniques are used between two staggered grids Paruta  _et al._ (2018). To
improve the numerical stability of our numerical simulations, parallel and
perpendicular numerical diffusions, such as

$\displaystyle
D(f)=\eta_{\perp}\left(\partial_{xx}^{2}+\partial_{yy}^{2}\right)f+\eta_{z}\partial_{zz}^{2}f,$
(24)

where $f$ denotes one of the evolved quantities, are added to the right-hand
side of all equations. We choose the perpendicular and parallel diffusion
coefficients, $\eta_{\perp}$ and $\eta_{z}$, to be constant and sufficiently
small not to affect significantly the results. The model is implemented in a
Fortran code using a MPI domain decomposition in all spatial directions. As a
consequence,, the computation cost of our simulations scales approximatively
linearly with the number of GMs taken into account.

The initial conditions of the turbulent nonlinear simulations impose equal
electron and ion densities and temperatures, such that
$n_{e}=\mathcal{N}^{00}$ and $T_{e}=T_{i}$ with top-hat-like profiles in the
perpendicular plane and uniform in $z$, similar to Eq. (5), such that
$f=\mathcal{A}_{0}[1-\tanh\left\\{(r-r_{s})/(4L_{s})\right\\}]/2$. We use
$\mathcal{A}_{0}=0.4$ for $f=n_{e}$ and $\mathcal{A}_{0}=0.4$ for $f=T_{e}$.
In addition, we set $\phi=\Lambda T_{e}$ to avoid unphysical and large
electron current into the sheath region. The initial values of
$\mathcal{N}^{20}$ and $\mathcal{N}^{01}$, given the initial ion density
$\mathcal{N}$ and ion temperature $T_{i}$ profiles, are obtained by inverting
Eq. (11), which yields $\mathcal{N}^{20}=N_{i}(T_{i}-1)/\sqrt{2}$ and
$\mathcal{N}^{01}=N_{i}(1-T_{i})$, with $T_{\parallel i}=T_{\perp i}=T_{i}$.
Finally, the parallel velocities, $U_{\parallel i}$ and $U_{\parallel e}$, are
initialized with smooth profiles along $z$, interpolating the values at the
end plates fixed according to the boundary conditions given in Eq. (23).
Random noise (with constant amplitude $0.01$) is added to the initial
densities, $n_{e}$ and $N_{i}$, temperatures $T_{i}$ (and associated GMs
$\mathcal{N}^{01}$, $\mathcal{N}^{02}$), and parallel velocities,
$U_{\parallel e}$ and $U_{\parallel i}$, profiles to seed turbulence.
Typically, a quasi-steady state is achieved after $100$ $c_{s0}/R$ time unit
(corresponding to $t\sim 4$ ms), where the sources of particle and energy are
compensated by the losses at the end plates.

##  IV Full-F Turbulent simulation results

In this section, we present the first turbulent and full-F simulations of the
GM approach of a linear plasma device, focusing on the parameters of the LAPD
experiment. We perform a comparison between the turbulent predictions of the
full-F GM approach (see Sec. II.2), with different numbers of GMs and values
of collisionality and compare them with the Braginskii model introduced in
Sec. II.6.

Our simulations parameters are similar to those used in Ref. Rogers and Ricci,
2010, where a helium LAPD plasma is considered. These parameters are sumarized
as follows: $n_{e0}=2\times 10^{12}$ cm-3, $T_{e0}=6$ eV, $T_{i0}=3$ eV
($\tau_{i}=0.5$), $\Omega_{i}\sim 960$ kHz, $\rho_{s0}=1.4$ cm,
$c_{s0}=1.3\times 10^{6}$ cm s-1, $m_{i}/m_{e}=400$ and $\nu_{0}=0.03$. The
LAPD vacuum chamber has a radius $R\simeq 0.56$ m (i.e., $R\simeq
40\rho_{s0}$) and a parallel length of $L_{z}\simeq 18$ m, such that we use
$L_{x}=L_{y}=100\rho_{s0}$ (or $L_{x}\sim L_{y}\sim 1.4$ m) and $L_{z}=36R$.
The reference time is $R/c_{s0}\sim 43$ $\mu$s. We consider the following
parameters for the density and temperature sources $L_{s}=1\rho_{s0}$,
$r_{s}=20\rho_{s0}$, $\mathcal{A}_{\mathcal{N}0}=\mathcal{A}_{T_{e0}}=0.04$
(with $\mathcal{A}_{\mathcal{N}\infty}=\mathcal{A}_{T_{e}\infty}=0.001$) and
$\mathcal{A}_{\mathcal{E}0}=0.02$ (with
$\mathcal{A}_{\mathcal{E}\infty}=\mathcal{A}_{\mathcal{N}\infty}$). We use a
numerical resolution of $N_{x}=N_{y}=192$ in the perpendicular plane, which
corresponds to a grid spacing of $\sim 0.52\rho_{s0}$. This resolution is also
sufficient to resolve turbulence occurring at perpendicular wavelengths longer
than $\rho_{s0}$ in the perpendicular plane, which is dominated by the long-
wavelength KH instability Rogers and Ricci (2010). On the other hand, a
coarser resolution of $N_{z}=64$ is used in the parallel direction. Given our
numerical resolution, we use $\eta_{\perp}=0.05$ and $\eta_{z}=5$ for
numerical stability. We remark that further numerical studies are required to
fully assess the impact of the numerical diffusion and resolution we use in
our simulations, but the values considered here are sufficient for the goals
of the present work.

In order to investigate the impact of ion collisions, we conduct a set of
nonlinear simulations in the high (HC) and low (LC) ion collisionality regime.
For each set, we consider different numbers of GMs $(P,J)$ to investigate the
convergence of the GM approach. More precisely, we consider $(P,J)=(2,1)$,
$(6,1)$, and $(12,1)$ in the LC regime and $(P,J)=(2,1)$ and $(6,1)$ in the HC
regime. We change the ion collisionality by varying the ion collision
frequency $\nu_{i}$ as an independent parameter while keeping all other
parameters constant. In the HC regime, the ion collision frequency is computed
using the LAPD physical parameters, such that
$\nu_{i}=1.38\sqrt{m_{i}/m_{e}}\nu_{0}/\tau_{i}^{3/2}\simeq 2.34$. In this
regime, the ion mean-free-path, $\lambda_{mpf}$, is considerably shorter than
the total length $L_{z}$, i.e.
$\lambda_{mpf}/L_{z}\simeq\sqrt{2\tau_{i}}R/L_{z}/\nu_{i}\ll 1$, and the
effects of the collision operator are expected to be important. On the other
hand, we set the ion collision frequency to be small in the LC regime, such
that $\nu_{i}\simeq 4\times 10^{-3}$ yielding $\lambda_{mpf}/L_{z}\sim 6.9$.
In this regime, the effect of the collision operator on the GMs is expected to
be negligible. In all cases, no artificial diffusion in velocity-space is
added, since the damping of the GMs by the Dougherty collision operator given
in Eq. (II.2) is sufficient to ensure stability of our simulations.

We remark that using $J=1$ is expected to be sufficient to represent the ion
distribution function $F_{i}$ in $x_{i}$ since kinetic effects such as
magnetic drifts (Frei, Hoffmann, and Ricci, 2022) and trapped particles (Frei
_et al._ , 2023) (which couple GMs with different $j$) are absent in the LAPD
configuration. In addition, neglecting the gyro-average in our ion model
allows us to remove the coupling between different $j$ driven by the
$x_{i}$-dependence of the gyro-average operator (Frei, Jorge, and Ricci,
2020). We also note that a number $J=1$ of Laguerre polynomials fully resolves
the density and energy sources (see Eq. (14)). As a matter of fact, the only
coupling to the $j\geq 2$ GMs is from the pitch-angle scattering contained in
the collision operator given in Eq. (II.2), which is dominated by collisional
damping. Indeed, if collisions are neglected GMs with $j>1$ are completely
decoupled from the the GMs with $j\leq 1$

This section is structured as follows. First, Sec. IV.1 provides an analysis
and comparison of simulations based on the full-F GM hierarchy, the cold-ion,
and the Braginskii models. Second, the turbulence characteristics are analyzed
and compared in more details in Sec. IV.2. Finally, we investigate the ion
distribution function in velocity-space in Sec. IV.3 and the GM spectrum in
quasi-steady state in Sec. IV.4 as a function of the number of GMs and for the
two collisionality regimes.

###  IV.1 Simulation results

This section presents a set of nonlinear and turbulent simulations of the LAPD
using the full-F GM hierarchy equation given in Eq. (II.2), the cold-ion model
in Eq. (17), and the Braginskii model introduced in Eq. (22).

A typical nonlinear evolution of the electron density, $n_{e}$, obtained by
using the GM hierarchy equation with $(P,J)=(6,1)$ GMs in the HC regime is
shown in Fig. 1. For $t\lesssim 21R/c_{s0}$, the profiles build up because of
the localized particle and energy sources present in the system. The steep
density and temperature gradients near $r\sim r_{s}$ drive an unstable
resistive drift-wave (not visible in Fig. 1 due to their small relative
amplitude), with the most unstable mode occurring at $k_{\perp}\rho_{s0}\sim
0.5$ ($k_{\perp}$ is the perpendicular wavenumber) with finite parallel
wavenumber and rotating in the ion diamagnetic direction. Large poloidal
flows, with velocity typically larger than the phase-velocity of the resistive
drift waves Popovich  _et al._ (2010), nonlinearly trigger a KH instability,
characterized by a long perpendicular wavelength ($k_{\perp}\rho_{s0}\ll 1$)
and $k_{\parallel}\simeq 0$. The KH instability becomes clearly visible around
$t\simeq 45R/c_{s0}$. This instability, which has been shown to dominate the
radial transport in LAPD Rogers and Ricci (2010), saturates at $t\simeq
53R/c_{s0}$, transporting the plasma to the $r\gtrsim r_{s}$ region and
yielding the broadening of the initial profiles. The role of the KH-dominated
transport in our simulations is confirmed by the strong steepening of the
profiles (similar to Ref. Rogers and Ricci, 2010) when the nonlinear term
$\left[\phi,\Omega\right]$ in Eq. (21) is artificially suppressed.

After $t\sim 130R/c_{s0}$ (similarly to previous GK and Braginskii turbulent
simulations of LAPD Popovich  _et al._ (2010); Pan  _et al._ (2018)), a quasi-
steady state is reached. At this point when the sources are compensated by the
losses at the end plates. In particular, the quasi-steady state in our
simulations can be monitored by examining the radial profiles. This is shown
in Fig. 2, where the initial radial profile of $n_{e}$ (associated with Fig.
1) reaches a quasi-steady state after $t\sim 130R/c_{s0}$. To support the
steady-state assumption in our simulations, we present the time-averaged
radial density profiles obtained over two time windows along with their
associated standard deviations. As observed, the time-averaged density
profiles remain the same (within their standard deviations). A similar
qualitative evolution is observed with a higher number of GMs in the LC and HC
regime, as well as in the cold-ion and Braginskii simulations.

Regarding the computational cost, we remark that reaching a quasi-steady state
requires a total of $\sim 41500$ CPU hours on the multi-core partition of the
Piz Daint supercomputer, approximatively.

Figure 1: Typical snapshots during the time evolution (from top-left to
bottom-right) of a LAPD simulation using $(P,J)=(6,1)$ GMs in the HC regime.
The normalized electron density $n_{e}$ (taken at $z=0$) is shown. First, the
top-hat-like sources fuel the central region ($r\lesssim r_{s}$) while a
resistive drift-wave instability is triggered as the gradients build up (top
left). Second, the KH instability develops (top right), saturates nonlinearly
(bottom left), and broadens the profiles to reach a quasi-steady state (bottom
right). $n_{e}$ is normalized, as well as all the quantities shown in the
following figures. Figure 2: Saturation of the radial profile of $n_{e}$
(binned by radius) associated with Fig. 1 at $z=0R$. Instantaneous radial
profiles are shown by the dashed lines for different times, and time-averaged
radial profiles over the time-windows $t\in[150,170]$ (red) and
$t\in[170,190]$ (blue) are also presented for comparison. The shaded colored
areas represent the standard deviations of each time-averaged profile. A
quasi-steady state is reached for $t\gtrsim 130$. The time $t$ is normalized
to $R/c_{s0}$.

The dynamics in the direction parallel to the magnetic field is shown in Fig.
3 during the quasi-steady state. Instantaneous snapshots of the parallel
turbulent structures of the electrostatic potential $\phi$, electron density
$n_{e}$, and electron temperature $T_{e}$ reveal elongated structures. All
quantities show weak parallel gradients $k_{\parallel}\simeq 0$, but have
slightly larger values at the center ($z=0$) and decrease in amplitude slowly
near the end plates (located at $z=-18R$ and at $z=18R$) due to the particle
and energy losses caused by the sheath boundary conditions. Similar parallel
structures are obtained when using a larger number of GMs. The turbulent
structures observed in Fig. 3 are in qualitative agreement with previous fluid
Rogers and Ricci (2010); Popovich  _et al._ (2010); Fisher  _et al._ (2015)
and GK Shi  _et al._ (2017); Pan  _et al._ (2018) turbulent simulations of
LAPD.

Figure 3: Typical snapshots of the parallel turbulent structures of the
electrostatic potential $\phi$ (top), electron density $n_{e}$ (middle), and
electron temperature $T_{e}$ (bottom) using $(P,J)=(6,1)$ GMs in the HC regime
during the quasi-steady state and as a function of $r$. We notice that
$\phi\sim 3T_{e}$.

We now examine the time-averaged radial profiles. These profiles are obtained
by averaging over a time window of $\sim 2$ ms during the quasi-steady state
as well as over the central region of LAPD, i.e. $-8R\leq z\leq+8R$ (or
$-4\text{m}\lesssim z\lesssim 4\text{m}$), a region commonly considered to
present experimental data Carter (2006) (a similar approach is used in
previous GK simulations Shi  _et al._ (2017); Pan  _et al._ (2018)). More
precisely, the averaged profile of the electrostatic potential $\phi$ is
computed by defining
$\left<\phi\right>_{t,z}=\int_{-\ell_{z}/2}^{+\ell_{z}/2}dz\int_{0}^{\tau}dt\phi/(\ell_{z}\tau)$,
with $\ell_{z}$ and $\tau$ are the $z$ and time length of the average windows,
respectively. The results are shown in Fig. 4, which displays the profiles of
$\left<\phi\right>_{t,z}$, $\left<n_{e}\right>_{t,z}$, and
$\left<T_{e}\right>_{t,z}$ obtained in the GM simulations, using different
numbers of GMs, in the cold-ion and the Braginskii simulations and as a
function of the radius. Instantaneous profiles are also included for
comparison.

We note, first, that the plasma profiles extend beyond $r_{s}$ illustrating
the broadening caused by the KH instability Rogers and Ricci (2010). More
precisely, the profiles are approximatively constant close to the center of
the device ($r<r_{s}$) and far from the source region ($r>r_{s}$), showing a
region of steep gradients near $r\sim r_{s}$, where the fluctuation level is
large (see Sec. IV.2). Second, the time-averaged radial profiles from the GM
simulations are very similar to the ones obtained from the Braginskii model.
Third, no noticeable differences are found between the simulations in the LC
and HC regimes and with different numbers of GMs. This suggests that ion
kinetic effects may not significantly influence the predictions of the
equilibrium (time-averaged) profiles in LAPD. On the other hand, the cold-ion
model consistently predicts larger time-averaged radial profiles, while the
gradients (not shown here) are of the same order as those obtained in the GM
and Braginskii simulations. This difference between the cold-ion model and
others cases with finite ion temperature is a consequence of the boundary
conditions given in Eq. (23), as discussed below. Fourth, the analysis of the
instantaneous profiles (indicated by dotted lines in Fig. 4) shows the
existence of large perpendicular turbulent structures associated with the KH
instability. Finally, we note that the time-averaged profiles obtained in Fig.
4 feature a similar behaviour of those obtained in previous fluid simulations
Rogers and Ricci (2010) and GK simulations Shi  _et al._ (2017); Pan  _et al._
(2018).

Figure 4: Time-averaged profiles (solid colored lines) of the electrostatic
potential $\left<\phi\right>_{t,z}$ (left), electron density
$\left<n_{e}\right>_{t,z}$ (center), and electron temperature
$\left<T_{e}\right>_{t,z}$ (right) obtained in the GM simulations in the HC
regime using $(P,J)=(2,1)$, $(6,1)$ and in the LC regime using $(P,J)=(2,1)$,
$(6,1)$ and $(12,1)$, as a function of the radial coordinate $r$. The cold-ion
and Braginskii simulations are also shown. The profiles are averaged over the
region $-8R\leq z\leq 8R$. The instantaneous profiles are shown by the dotted
colored lines for the Braginskii and $(P,J)=(6,1)$ GM in the HC regime
simulations. The gray shaded areas represent the radial extent of the
localized sources.

We remark that the electrostatic potential profile $\phi$ follows
approximatively the electron temperature $T_{e}$, as shown in Fig. 4. Indeed,
$\phi\sim\Lambda T_{e}$ is required to have comparable electron and ion
outflows in steady-state, such that $U_{\parallel i}\sim U_{\parallel e}$ near
the end plates, according to Eq. (23). To verify that $\phi\sim\Lambda T_{e}$
in our simulations, we evaluate the radial profile of the instantaneous
difference, $\phi-\Lambda T_{e}$, taken at the center of the device ($z=0R$)
for the GM (both the LC and HC regimes are considered), cold-ion and
Braginskii simulations during the quasi-steady state. The results are shown in
Fig. 5. We first observe that the GM and Braginskii simulations yield similar
negative $\phi-\Lambda T_{e}$ values. On the other hand, $\phi-\Lambda T_{e}$
is roughly constant and approximatively vanishes for all radii in the cold-ion
model. The deviations observed in $\phi-\Lambda T_{e}$ are the result of the
boundary conditions specified in Eq. (23). In fact, to achieve an equilibrium
between parallel ion and electron flows, it is necessary that $\phi-\Lambda
T_{e}\simeq-T_{e,s}\ln\left(1+\tau_{i}T_{i,s}/T_{e,s}\right)/2$ (as indicated
by the black dotted line in Fig. 5) when finite ion temperature is considered
in the case of Braginskii and GM simulations. In contrast, this condition
reduces to $\phi-\Lambda T_{e}\simeq 0$ in the cold-ion model. Consequently,
due to the finite ion temperature dependence of the boundary condition Eq.
(23b), the equilibrium electrostatic potential $\phi$ exhibits higher
amplitude in the cold ion model, which is confirmed in the time-averaged
profiles observed in Eq. (4). Finally, in all cases considered, $\phi-\Lambda
T_{e}$ remains smaller than the values of $\phi$ and $\Lambda T_{e}$
($\phi-\Lambda T_{e}\sim 0.1$ for $r\lesssim r_{s}$ compared to
$\phi\sim\Lambda T_{e}\sim 2$, see Fig. 3).

Figure 5: Radial instantaneous profiles of the difference between $\phi$ and
$\Lambda T_{e}$ at $z=0R$, in the LC and HC regimes for different number of
GMs $(P,J)$ and at quasi-steady state. The results from the cold-ion and
Braginskii simulations are also shown. The black dotted line shows the value
of $-T_{e,s}\ln\left(1+~\tau_{i}T_{i,s}/T_{e,s}\right)/2$ in the Braginskii
simulation at $z=0R$ (similar values of obtained in the GM simulations). The
gray shaded area ($r\leq r_{s}$) represents the radial extent of the top-hat-
like sources. For $r\gg r_{s}$, $\phi$ relaxes to $\Lambda T_{e}$.

###  IV.2 Turbulence analysis

Figure 6: Snapshots of the electrostatic potential $\phi$ (top) and of the
fluctuations $\widetilde{\phi}=\phi-\left<\phi\right>_{t}$ (bottom) taken at
$z=0R$ obtained in the Braginskii, cold-ion (cold-$T_{i}$), $(P,J)=(2,1)$,
$(12,1)$ in the LC regime, and $(P,J)=(6,1)$ in the HC regime simulations
(from left to right).

We now delve into the analysis of the turbulence properties, comparing the GM
predictions with the Braginskii simulations. The instantaneous fluctuations
are obtained by subtracting the time-averaged profiles from the full
quantities, such that the fluctuation of, e.g., the electrostatic potential,
$\widetilde{\phi}$, is defined by
$\widetilde{\phi}=\phi-\left<\phi\right>_{t}$, where
$\left<\phi\right>_{t}=\int_{0}^{\tau}dt\phi/\tau$ denotes the time-averaged
potential. Similar definitions for the other quantities are used. The top
panels of Fig. 6 show instantaneous snapshots of $\phi$ in the plane
perpendicular at the center of the device $z=0R$, while the bottom panels
illustrate $\widetilde{\phi}$ snapshots. The Braginskii, cold-ion, and GM
simulations with various $(P,J)$ are considered.

Figure 7: Snapshots of the electron density $n_{e}$ (top) and of the
fluctuation $\widetilde{n_{e}}=n_{e}-\left<n_{e}\right>_{t}$ (bottom), similar
to Fig. 6.

We first observe that the fluctuations in the Braginskii model exhibits
similar structures than those obtained in Ref. Fisher  _et al._ , 2015. In
particular, the level of fluctuations is low at the center of the device and
far from the source region, while it is large where the equilibrium gradient
is steeper, in particular near $r\sim r_{s}$ (see Fig. 9). Notably, the
$\widetilde{\phi}$ snapshots reveal the presence of large amplitude structures
propagating outwards. These observations hold for all the GM simulations,
demonstrating qualitative agreement between the GM approach and the Braginskii
model. While the fluctuations of the potential $\phi$ is not significantly
affected by the number of GMs used in the simulation or by the collisionality
regime, pointing out the fact that the KH instability (which drives
turbulence) has a fluid nature, minor differences in the turbulent properties
can still be observed. In fact, the use of a small number of GMs tends to
produce slightly larger turbulent structures. This can be observed, for
instance, by comparing the results of $(P,J)=(2,1)$ with the $(12,1)$
simulations in the LC regime. Finally, we observe that the cold-ion model
produces the largest turbulent structures, which is consistent with the broad
time-averaged profiles observed in Fig. 4. The same observations apply to the
snapshots of the electron density $n_{e}$ and its associated fluctuations
$\widetilde{n_{e}}$, as shown in Fig. 7. Similar plots are obtained for
$T_{i}$ and $T_{e}$, but not shown. Finally, it is worth noting that the
fluctuations of the electrostatic potential $\widetilde{\phi}$ normalized to
the electron temperature $T_{e}$ can be of the order of unity, as inferred
from Eq. (4) and Eq. (6). In particular, we observe
$\widetilde{\phi}/T_{e}\lesssim 1$ at all radii with, for example,
$\widetilde{\phi}\sim 0.5$ and $T_{e}\sim 0.6$ for $r\lesssim r_{s}$.

In Fig. 8, we present snapshots of the vorticity for the $\bm{E}\times\bm{B}$
flow, which is described by the vorticity variable
$\Omega=\gradient_{\perp}^{2}\phi$, defined in Eq. (21), at the center of the
device ($z=0R$) at steady-state obtained in the case of the GM simulations
with $(P,J)=(6,1)$ in the LC and HC regimes. It is observed that the vorticity
exhibits large perpendicular structures, i.e. $k_{\perp}\rho_{s0}\ll 1$,
driven by the KH instability, which develop in the region of steep equilibrium
gradients near $r\sim r_{s}$. These structures have slightly smaller amplitude
away from $r_{s}$ in all cases. Additionally, no clear differences are
observed between the LC and HC regimes. Similar plots are obtained with the
Braginskii and GM simulations. Finally, we note that further numerical
investigations are necessary to verify the numerical resolution of the
smallest perpendicular scale revealed by Fig. 8.

Figure 8: Snapshots of the vorticity variable, $\Omega$, obtained in the GM
simulations with $(P,J)=(6,1)$ in the LC (left) and HC (right) regimes at
$z=0R$. Figure 9: RMS of the electron density $n_{e}$ (normalized to
$\left<n_{e}\right>_{t}(r)$) (top) and electrostatic potential $\phi$
(normalized to $\left<\phi\right>_{t}(r)$) (bottom) fluctuations (computed
from Fig. 7 and Fig. 6) as a function of the radius $r$ in the case of the
Braginskii, cold-ion, and GM simulations at $z=0R$. The shaded area indicates
the radial extent of the sources.

We now proceed to analyze the root mean square (RMS) of the fluctuations,
defined as $\sqrt{\left<\widetilde{n_{e}}^{2}\right>_{t}}$ in the case of the
electron density fluctuations $\widetilde{n_{e}}$ (and similarly for the other
quantities). Fig. 9 displays the RMS of the electron density, $n_{e}$, and the
electrostatic potential, $\phi$, fluctuations plotted as a function of the
radius. The data are computed at $z=0R$ and normalized to
$\left<n_{e}\right>_{t}(r)$ (and to $\left<\phi\right>_{t}(r)$) Friedman  _et
al._ (2012); Shi  _et al._ (2017). We find that the RMS values of the density
displayed in Fig. 9 exhibit a qualitatively similar behaviour to those
obtained in previous fluid simulations Fisher  _et al._ (2015); Ross  _et al._
(2019) and GK simulations Shi  _et al._ (2017). Consistent with the
observations made in Fig. 7, the RMS values reach their maximum when the
gradients are most pronounced near $r\sim r_{s}$. For $r\lesssim r_{s}$ and
for $r\gtrsim r_{s}$ (where the gradients are smaller), the RMS values
decrease because of the absence of the instability drive. Using a low number
of GMs or considering the LC regime results in slightly larger RMS values (in
particular of $\widetilde{\phi}$). Overall, this indicates that the level of
fluctuations in the steep gradient region can be sensitive to the number of
GMs used in the simulations. Finally, we note that the GM simulation with
$(P,J)=(6,1)$ in the HC regime is the closest to the Braginskii predictions,
and the largest RMS values (especially for $\widetilde{n_{e}}$) are obtained
in the cold-ion model.

We compare the RMS of the parallel electrical current $J_{\parallel}$ measured
at the sheath entrance located at $z=-18R$. The results are shown in Fig. 10
as a function of the radius and normalized to the maximum of
$\left<J_{\parallel}\right>_{t}(r)$. It is clearly observed that the boundary
conditions imposed on the electron and ion parallel velocities allow for the
parallel current to fluctuate. This is in contrast to the case of logical
sheath boundary condition, where $J_{\parallel}=0$ is imposed everywhere
Parker  _et al._ (1993). We remark that larger fluctuations of $J_{\parallel}$
are obtained in the Braginskii simulations, while the largest RMS is observed
in the case of the cold-ion model, a consequence of the finite ion temperature
effect in the boundary conditions retained in the former (see Eq. (23)) and
ignored in the later.

Figure 10: RMS of the total parallel current $J_{\parallel}$ (normalized to
the maximum of $\left<J_{\parallel}\right>_{t}(r)$) as a function of the
radius at the sheath entrance near $z=-18R$.

We now turn our attention to the skewness of the ion density fluctuations,
which is defined as the third normalized moment of the electron density
fluctuation, that is
$\left<\widetilde{n_{e}}^{3}\right>_{t}/\left<\widetilde{n_{e}}^{2}\right>_{t}^{3/2}$.
The skewness of the density is often used to characterize the dominance of
plasma holes and blobs, associated with negative and positive skewness
respectively Friedman  _et al._ (2012); Fisher  _et al._ (2015); Ross  _et
al._ (2019). Fig. 11 shows the skewness of the electron density $n_{e}$. In
all cases, the skewness is negative for $r\lesssim r_{s}$, indicating the
presence of density holes in the region where the plasma source is present. On
the other hand, in the region where $r\gtrsim r_{s}$, the skewness is
positive. The sign and amplitude of the skewness shown in Fig. 11 have similar
qualitative behaviour as in previous fluid Ross  _et al._ (2019) and GK Shi
_et al._ (2017); Pan  _et al._ (2018) simulations. In particular, the values
obtained in the GM simulations are of the same order to those observed in the
Braginskii case, albeit slightly smaller. Overall, the present turbulent
analysis demonstrates that the full-F GM approach is in qualitative agreement
with the Braginskii model, which is also employed in previous numerical
investigations Fisher  _et al._ (2015); Ross  _et al._ (2019) and validated
with experimental data Friedman  _et al._ (2012). We note that a more
comprehensive analysis of turbulence, including aspects such as particle and
energy transport, is necessary to evaluate the impact of velocity-space
resolution and collisions (along with collision operator models) in greater
detail. This analysis is deferred to future investigations, which would
involve extending the present work to include FLR effects and removing the
Boussinesq approximation, important for accurate turbulence transport
predictions.

Figure 11: Skewness of the electron density $n_{e}$ as a function of $r$ at
$z=0R$. The gray shaded area represents the radial extent of the sources.

###  IV.3 Ion distribution function at quasi-steady state

Figure 12: Quasi-steady state ion distribution function, $F_{i}$, plotted as a
function of $(v_{\parallel}/v_{Ti},x_{i})$, obtained using $(P,J)=(6,1)$ GMs
in the HC regime (top) and LC regime (bottom). $F_{i}$ is computed at the
$z=-18R$ (left), $z=0R$ (center), and $z=18R$ (right) at the center in the
perpendicular plane ($x=y=0$). Figure 13: Cuts of the ion distribution
function $F_{i}$ (normalized to its maximum) at $x_{i}=0$ at the right sheath
entrance near $z=18R$ and $x=y=0$ (top) for different $(P,J)$ in the HC (solid
lines) and LC (dotted lines) regimes, and close view over the region
$v_{\parallel}/v_{Ti}\in\left[1,2.5\right]$ (bottom). The solid black line
represents a Maxwellian distribution function shifted by the ion sound speed
$c_{s}$ (see Eq. (6)). A similar plot is obtained at the left sheath entrance
($z=-18R$).

We now investigate the features of the ion distribution function $F_{i}$ in
velocity-space. To obtain the full-F ion distribution function, $F_{i}$, from
the GM simulations, we use the expansion in Eq. (8), truncated to a finite
number of GMs, and we compute it as a function of $x_{i}$ and the unshifted
parallel coordinate, given by $v_{\parallel}/v_{Ti}=s_{\parallel
i}+\sqrt{2\tau_{i}}U_{\parallel i}$. Also for this analysis, we consider the
quasi-steady period. Figure. 12 shows $F_{i}$ obtained from the $(P,J)=(6,1)$
simulations in the HC and LC regime at the center of the machine ($z=0R$) and
at the sheath entrances, $z=-18R$ and $z=18R$. At the two sheath entrances,
$F_{i}$ is centered around the ion parallel velocity, $U_{\ i}=\pm c_{s}$
respectively, a consequence of the Bohm sheath boundary conditions given in
Eq. (23). On the other hand, $F_{i}$ is centered around $v_{\parallel}\simeq
0$ at $z=0R$, where $U_{\parallel i}\simeq 0$. The absence of fine velocity-
space structures in Fig. 12 in both the HC and LC regimes confirms the weak
dependence of turbulence properties on the number of GMs, as reported in Sec.
IV.2. In fact, fine velocity-space structures in $F_{i}$ are not expected to
play a significant role in LAPD due to the absence of magnetic drifts, trapped
particles, of the high-collisional regime, and the fluid nature of the KH-
dominated turbulent transport. Nevertheless, we remark that strong velocity-
space gradients can appear near the sheath entrance in the electron
distribution function (see, e.g., Ref. Pan  _et al._ , 2018). The Maxwellian
shape of $F_{i}$ observed in our simulations is also consistent with the fact
that both two-fluid Braginskii Rogers and Ricci (2010); Ross  _et al._ (2019)
and GK Shi  _et al._ (2017); Pan  _et al._ (2018) simulations of LAPD yield
qualitatively similar results.

Figure. 13 shows the ion distribution function at the sheath entrance ($z=18R$
and $x=y=0$) for $x_{i}=0$, in the LC and HC regimes and for different values
of $(P,J)$. We first observe that the bulk region of $F_{i}$ (near
$v_{\parallel}/v_{Ti}\sim 1$) is well approximated by a shifted Maxwellian.
However, deviations from the Maxwellian distribution function are noticeable
in the tails of $F_{i}$ in the LC regime. These deviations become pronounced
as $(P,J)$ increases (e.g., from $(6,1)$ to $(12,1)$) which indicates that
$F_{i}$ is not sufficiently resolved in the LC regime. Finally, we remark that
collisional effects tend to widen $F_{i}$ due to the collisional parallel
velocity-space diffusion present in the nonlinear Dougherty operator. We note
that further convergence studies are required to fully evaluate the velocity-
space resolution in Fig. 13 by considering a larger number of GMs than those
considered in the present work, in particular in the LC regime.

The GM hierarchy equation in Eq. (II.2) does not enforce the positivity of the
ion distribution function when evolving a finite number of GMs. For instance,
the use of $v_{\parallel}/v_{Ti}$ as an argument in the Hermite polynomials,
$H_{p}$, in Eq. (8) would compromise the convergence properties of the GM
approach, with respect to the use of $v_{\parallel}/v_{Ti}$, leading to
simulations that show unphysical distribution functions with negative values
when the same number of GMs are considered, as in the simulations presented
here (Fig. 12 would show negative values). In fact, if the unshifted GMs
$\mathcal{N}_{v_{\parallel}}^{pj}$, defined with respect to
$v_{\parallel}/v_{Ti}$ as the argument of $H_{p}$, i.e.

$\displaystyle\mathcal{N}^{pj}_{v_{\parallel}}=2\pi\int_{-\infty}^{\infty}dv_{\parallel}\int_{0}^{\infty}d\mu\frac{B}{m_{i}}F_{i}\frac{H_{p}\left(\frac{v_{\parallel}}{v_{Ti}}\right)L_{j}(x_{i})}{\sqrt{2^{p}p!}},$
(25)

are used to expand $F_{i}$, it is found that
$\mathcal{N}_{v_{\parallel}}^{pj}\neq 0$ for $(p,j)>0$, even when $F_{i}$ is a
Maxwellian distribution function centered at $U_{\parallel i}\neq 0$. Indeed,
using Eq. (25), one derives the analytical expression of the unshifted GMs for
$F_{i}=F_{Mi}$,

$\displaystyle\mathcal{N}_{v_{\parallel}}^{pj}$
$\displaystyle=\frac{\delta_{j}^{0}}{\sqrt{\pi}}\int_{-\infty}^{\infty}d\left(\frac{v_{\parallel}}{v_{Ti}}\right)e^{-(v_{\parallel}/v_{Ti}-U_{\parallel
i}/(\sqrt{2\tau_{i}}))^{2}}\frac{H_{p}\left(\frac{v_{\parallel}}{v_{Ti}}\right)}{\sqrt{2^{p}p!}}$
$\displaystyle=\sqrt{\frac{2^{p}}{p!}}\left(\frac{U_{\parallel
i}}{\sqrt{2\tau_{i}}}\right)^{p}\delta_{j}^{0},$ (26)

where $U_{\parallel i}$ is normalized to $c_{s0}$. While the amplitude of the
unshifted GM decreases rapidly in the presence of subsonic ion flow,
$U_{\parallel i}\ll 1$, the decrease of the amplitude with $p$ is slower in
the presence of sonic flows, such that
$\mathcal{N}^{pj}_{v_{\parallel}}\sim\sqrt{2^{p}/p!}$.

###  IV.4 GM spectrum at quasi-steady sate

Figure 14: Absolute value of the GMs $\mathcal{N}^{pj}$ (normalized to
$\mathcal{N}^{00}$) associated with the distribution functions shown in Fig.
13, plotted on a logarithmic scale as a function of $p\geq 2$ (solid lines for
$j=0$ and dashed lines for $j=1$) obtained from different $(P,J)$ simulations
in the LC (square colored markers) and HC (circle colored markers) regime at
$z=18R$ and $x=y=0$. The amplitude of the GM decreases with $p$ in all cases
and it is faster in the HC regime.

To better assess the velocity-space representation of $F_{i}$ in our
simulations, we plot the absolute value of the GMs, $\mathcal{N}^{pj}$
(normalized to $\mathcal{N}^{00}$), at the sheath entrance of the device,
$z=18R$ and $r=0$, in Fig. 14. We consider the GMs associated with the
distribution functions displayed in Fig. 13. As it can be clearly observed,
the amplitude of the GMs decays faster in the HC regime than in the LC regime.
In fact, extrapolating from linear theory in slab geometry Frei, Hoffmann, and
Ricci (2022), the amplitude of the GMs is expected to follow
$\left|\mathcal{N}^{pj}\right|\sim\exp\left[-c_{p}(\nu_{i})p^{\alpha}\right]$
in the $p\gg 1$ limit (with $c_{p}(\nu_{i})$ a positive coefficient that
increases with collisionality, $\nu_{i}$ and $\alpha>0$). While this estimate
has been verified in linear studies and remains to be proven valid in full-F
turbulent simulations, it provides us with a useful proxy for the expected
decay of the GMs in the presence of collisions.

The results show that $P>12$ would be required to ensure that $F_{i}$ is
sufficiently well resolved in the LC regime. Since the smaller relevance of
the LC regime for the experimental parameters of LAPD, simulations with $P>12$
have not been carried out in the present work to properly assess the velocity-
space resolution of $F_{i}$. On the other hand, the contributions from
$\mathcal{N}^{p0}$ with $p\gtrsim 4$ are negligible in the HC regime, thereby
justifying the closure by truncation for $P\gtrsim 4$. We also notice that
$\mathcal{N}^{10}=0$ in all cases, as a consequence of Eq. (11b). Finally, we
note that the amplitude of the low-order GMs is not sensitive to $P$, as shown
in Fig. 14. More precisely, the low-order GMs for $(P,J)=(6,1)$ strongly
resemble the ones of the $(P,J)=(12,1)$ simulation in the LC case. This holds
true also in the HC regime, for instance, by comparing the $(P,J)=(6,1)$ and
$(P,J)=(2,1)$ simulations.

We also show in Fig. 14 that the amplitudes of the $j=1$ GMs are smaller than
the ones with $j=0$. Numerical experiments (not shown) have also revealed that
GMs with $j>1$ have a negligible amplitude due to the absence of source for
$j>1$ and of FLR effects. This suggests (in addition to the similar results
obtained in Sec. IV.2 with different $(P,J)$) that full-F turbulent
calculations using the GM approach are less sensitive to the values of $P$ and
$J$ than linear computations Jorge  _et al._ (2019); Frei, Hoffmann, and Ricci
(2022), where applying a closure by truncation at low $P$ and $J$ can
introduce spurious artifacts Frei  _et al._ (2023). Otherwise, Fig. 14 reveals
that the slightly larger RMS values depicted in Fig. 9 (e.g., $(P,J)=(6,1)$ in
the LC and $(P,J)=(2,1)$ in the HC regime) correspond to cases where the GM
representation of $F_{i}$ is unresolved. Additional investigations are
required to verify the effect of closure in the presence of kinetic effects
such as trapped particles and magnetic drifts, which are absent in LAPD, and
of FLR effects, which are neglected in our model.

Finally, Fig. 15 presents snapshots of the GMs for different values of $p$ in
the perpendicular plane obtained for the $(P,J)=(6,1)$ simulations in the LC
and HC regimes. It is also visible that the turbulent structures are dominated
by a long-wavelength perpendicular fluctuations for all values of $p$, which
are driven by the KH instability. Indeed, the amplitude of these turbulent
structures is strongly reduced if the nonlinear drive of the KH instability is
artificially removed from our simulations. The decay of the amplitude of the
turbulent structures due to collisions and with increasing $p$ is also
evident.

Figure 15: Instantaneous snapshots of the GMs $\mathcal{N}^{pj}$, with
$p=2,3,4,6$ and $j=0$, in the $(x,y)$ perpendicular plane at the sheath
entrance located at $z=18R$ from the $(P,J)=(6,1)$ simulations in the LC (top)
and HC (bottom) regime. The large perpendicular structures are caused by the
KH instability.

##  V Conclusions

In this work, we present the first full-F turbulent simulations based on the
GM approach using an arbitrary number of moments in a linear plasma device
configuration with open straight field lines, such as LAPD. Neglecting FLR
effects, we consider an electrostatic and long-wavelength ion GK model for the
full ion distribution function $F_{i}$, coupled to the electron Braginskii
fluid model for the electron density $n_{e}$, parallel velocity $U_{\parallel
e}$, and temperature $T_{e}$. The ion GK model is solved by deriving a full-F
ion GM hierarchy equation, based on the Hermite-Laguerre polynomial expansion
of $F_{i}$. In particular, a velocity-space coordinate centered at the local
fluid ion parallel velocity is used to expand $F_{i}$, which ensures good
convergence properties of the Hermite expansion in the presence of sonic ion
flows. The GM hierarchy equation we consider is equivalent to the
electrostatic and Boussinesq limits of the GK moment model for the boundary
region derived in Ref. Frei, Jorge, and Ricci, 2020 with vanishing FLR
effects. To account for the parallel losses at the end plates, Bohm sheath
boundary conditions are used, equivalent to those previously used in LAPD
fluid simulations Rogers and Ricci (2010). A nonlinear ion-ion Dougherty
collision operator is also considered. The ion GM hierarchy equation with the
fluid electron model is implemented in a numerical code, allowing us to
perform the first full-F turbulent calculations based on a moment approach
with a flexible number of GMs. This is in contrast to previous full-F
gyrofluid simulations where the number of moments is fixed Wiesenberger  _et
al._ (2019); Galassi  _et al._ (2022); Wiesenberger and Held (2022).

We present the simulations of a linear device using LAPD physical parameters
based on a Helium plasma Rogers and Ricci (2010) and a first-of-the-kind
comparison with the two-fluid Braginskii model. Several nonlinear simulations
are performed using a different number of Hermite and Laguerre GMs in a low
and high-collisional ion regime. Overall, a qualitative agreement on the time-
averaged radial profiles between the Braginskii model and the GM approach is
observed. This is expected from our analysis which shows that turbulence is
dominated by the long perpendicular wavelength and $k_{\parallel}\simeq 0$
Kelvin-Helmoltz instability of fluid nature.

The RMS and skewness of the fluctuations in the GM simulations also agree with
the ones previously obtained in fluid Rogers and Ricci (2010); Fisher  _et
al._ (2015) and GK Shi  _et al._ (2017); Pan  _et al._ (2018) simulations of
LAPD. In particular, we find that the RMS values are often larger than the
ones predicted by the Braginskii model, if the number of GMs is not sufficient
to properly resolve the ion distribution function. The largest RMS values are
observed with the cold-ion reduced model (with a difference up to $\sim 20\%$
with respect to the Braginskii model), while the results closest to the one of
the Braginskii model are obtained if collisions are introduced in the GM
approach with a sufficient number of GMs, in this case $(P,J)=(6,1)$. Overall,
collisions reduce slightly the turbulent fluctuations level, but they do not
significantly alter the observed turbulent regimes. At the same time, the
analysis of the ion distribution function $F_{i}$ reveals that collisions damp
the amplitudes of the GMs, thereby allowing for a reduction in the number of
GMs required in the simulations (typically from $(P,J)\sim(12,1)$ in the low
collisional regime to $(P,J)\sim(6,1)$ in the high-collisional regime of
LAPD).

We note that although the present work allows the simulation of LAPD with a
flexible number of GMs, unlike previous full-F gyrofluid
simulationsWiesenberger  _et al._ (2019); Galassi  _et al._ (2022);
Wiesenberger and Held (2022), improvements in our physical model are required
to increase the fidelity of our simulations. For example, the inclusion of FLR
effects Held, Wiesenberger, and Kendl (2020); Held and Wiesenberger (2022) is
necessary to study small scale fluctuations. A more detailed study of the role
of FLR effects is the subject of future work. Incorporating FLR effects with a
flexible number of GMs necessitates conducting comprehensive numerical
investigations to evaluate the impact of the numerical and velocity-space
resolution on turbulent structures, in particular at small scales, supported
by the examination of conservation laws, an aspect not explored in the present
work. In addition, a more accurate description of the role of collisions
requires the implementation of a nonlinear collision operator model with
increasing physical fidelity, such as the nonlinear Coulomb operatorJorge,
Frei, and Ricci (2019). Proper sheath boundary conditions for the GM hierarchy
equation, which extend the simplified Bohm sheath boundary condition used here
(see Eq. (23)), can improve the reliability of our simulations. These boundary
conditions can be obtained by a procedure similar to that described in Ref.
Loizu  _et al._ , 2012. Finally, we note that the implementation of a kinetic
electron description is also essential to perform high-fidelity LAPD
simulations, since fast and less collisional electrons (with $T_{e}\sim 15$
eV) are emitted by pulsed plasma discharges in experiments Gekelman  _et al._
(2016). Furthermore, kinetic electrons are important in setting the sheath
boundary conditions where electrons are reflected due to the potential drop,
giving rise to strong speed-space gradients in the electron distribution
function Pan  _et al._ (2018).

More generally, the present work represents a step towards the development of
future full-F turbulent simulations of the boundary region of fusion devices
using the GM approach with a flexible number of GMs, which may provide an
ideal flexible tool to capture kinetic and collisional effects at the desired
level of accuracy.

## Acknowledgement

The authors acknowledge helpful discussions with Alessandro Geraldini and
Stephan Brunner. This work has been carried out within the framework of the
EUROfusion Consortium, via the Euratom Research and Training Programme (Grant
Agreement No 101052200 — EUROfusion) and funded by the Swiss State Secretariat
for Education, Research and Innovation (SERI). Views and opinions expressed
are however those of the author(s) only and do not necessarily reflect those
of the European Union, the European Commission, or SERI. Neither the European
Union nor the European Commission nor SERI can be held responsible for them.
The simulations presented herein were carried out in part on the CINECA
Marconi supercomputer under the TSVVT422 project and in part at CSCS (Swiss
National Supercomputing Center). This work was supported in part by the Swiss
National Science Foundation.

##  Appendix A High-collisional limit of ion temperature equation

In this appendix, we derive the Braginskii equation for the ion temperature
$T_{i}$, given Eq. (22b), from the ion GM hierarchy equation provided in Eq.
(II.2). This allows us to elucidate the relationship between the GM approach
and the Braginskii model, particularly in the high-collisional regime, and to
justify the terms on the right-hand side of Eq. (22b).

First, we note that since our model neglects FLR effects, the gyrocenter
quantities are assumed to be equal to the particle fluid quantities. Hence,
the time evolution of the ion temperature $T_{i}$ can be expressed in terms of
GMs using the definitions of the gyrocenter fluid quantities in Eq. (11).
Evaluating the time derivative of $T_{i}$ yields

$\displaystyle
N_{i}\partial_{t}T_{i}=\frac{\sqrt{2}}{3}\partial_{t}\mathcal{N}^{20}-\frac{2}{3}\partial_{t}\mathcal{N}^{01}+(1-T_{i})\partial_{t}N_{i}.$
(27)

The time derivatives of the GMs contained in Eq. (27) can be evaluated using
the ion GM hierarchy equation given in Eq. (II.2) with $(p,j)=(2,0)$ and
$(0,1)$. We obtain

$\displaystyle\partial_{t}\mathcal{N}^{20}$
$\displaystyle=-\sqrt{\tau_{i}}\partial_{z}\sqrt{3}\mathcal{N}^{30}-\partial_{z}\left(U_{\parallel
i}\mathcal{N}^{20}\right)-\frac{1}{\rho_{*}}\left[\phi,\mathcal{N}^{20}\right]$
$\displaystyle-\left(2\mathcal{N}^{20}+\sqrt{2}\mathcal{N}^{00}\right)\partial_{z}U_{\parallel
i}+\mathcal{C}^{20}_{i}+S_{\mathcal{E}}^{20},$ (28a)
$\displaystyle\partial_{t}\mathcal{N}^{01}$
$\displaystyle=-\sqrt{\tau_{i}}\partial_{z}\mathcal{N}^{11}-\partial_{z}\left(U_{\parallel
i}\mathcal{N}^{01}\right)-\frac{1}{\rho_{*}}\left[\phi,\mathcal{N}^{01}\right]$
$\displaystyle+\mathcal{C}^{01}_{i}+S_{\mathcal{E}}^{01}.$ (28b)

Using Eq. (28) with the ion continuity equation given in Eq. (17a) into Eq.
(27), we derive

$\displaystyle N_{i}\partial_{t}T_{i}$
$\displaystyle=-\frac{\sqrt{2\tau_{i}}}{3}\partial_{z}\left(\sqrt{3}\mathcal{N}^{30}-\sqrt{2}\mathcal{N}^{11}\right)$
$\displaystyle-U_{\parallel
i}N_{i}\partial_{z}T_{i}-\frac{N_{i}}{\rho^{*}}\left[\phi,T_{i}\right]-\frac{2}{3}P_{\parallel
i}\partial_{z}U_{\parallel i}$
$\displaystyle+S_{\mathcal{E}}+\left(1-T_{i}\right)S_{\mathcal{N}},$ (29)

where the collision terms, $\mathcal{C}^{20}_{i}$ and $\mathcal{C}^{01}_{i}$,
cancel out due to the energy conservation of the collision operator (see Eq.
(II.2)). In addition, the definitions of the sources given in Eq. (14) are
used to derive the last terms on the right-hand side of Eq. (A).

We now use the fact that FLR effect are neglected in our model such that
$N_{i}\simeq n_{e}$. Moreover, an expression for the parallel gradient of the
parallel ion flow, $\partial_{z}U_{\parallel i}$, can be obtained from the
vorticity equation given in Eq. (22). At the leading-order in the long-
wavelength limit, one finds that $\partial_{z}J_{\parallel}\simeq 0$, such
that

$\displaystyle N_{i}\partial_{z}U_{\parallel i}\simeq\left(U_{\parallel
e}-U_{\parallel i}\right)\partial_{z}n_{e}+n_{e}\partial_{z}U_{\parallel e}.$
(30)

Finally, introducing the normalized parallel ion heat flux,

$\displaystyle q_{\parallel
i}=\sqrt{2\tau_{i}}\left(\frac{\sqrt{3}}{2}\mathcal{N}^{30}-\frac{\mathcal{N}^{11}}{\sqrt{2}}\right),$
(31)

and using Eq. (30), the ion temperature equation Eq. (A) becomes

$\displaystyle\partial_{t}T_{i}+\frac{1}{\rho^{*}}\left[\phi,T_{i}\right]+U_{\parallel
i}\partial_{z}T_{i}$ $\displaystyle=+\frac{2}{3}T_{\parallel
i}\left[\left(U_{\parallel i}-U_{\parallel
e}\right)\frac{\partial_{z}n_{e}}{n_{e}}-\partial_{z}U_{\parallel e}\right]$
$\displaystyle-\frac{2}{3}\frac{1}{n_{e}}\partial_{z}q_{\parallel
i}+\frac{A_{\mathcal{E}}}{n_{e}}+\left(1-T_{i}\right)\frac{S_{\mathcal{N}}}{n_{e}}.$
(32)

We remark that Eq. (A) can also be deduced from other full-F gyrofluid models
Held  _et al._ (2016); Wiesenberger and Held (2022), which consider finite ion
temperature effects with a fixed number of GMs. To close Eq. (A), an
expression for the parallel ion heat flux $q_{\parallel i}$ must be
prescribed. In contrast to the ion GM hierarchy equation, which treats the ion
parallel heat flux $q_{\parallel i}$ as a dynamical variable by evolving high-
order GMs (see Eq. (31)), the Braginskii model assumes a high-collisional
limit (where $\lambda_{mpf}k_{\parallel}\ll 1$), allowing for the derivation
of a collisional closure for $q_{\parallel i}$. Specifically, using the
Fokker-Planck collision operator and applying the high-collisional closure via
the Chapman-Enskog procedure Chapman and Cowling (1941); Braginskii (1965);
Jorge, Ricci, and Loureiro (2017), we can obtain the fluid closure
$q_{\parallel i}\simeq-\chi_{\parallel}\partial_{z}T_{i}$, where
$\chi_{\parallel}$ represents the normalized parallel ion thermal conductivity
predicted by the Braginskii transport equation Braginskii (1965); Zeiler,
Drake, and Rogers (1997). Hence, in the high-collisional limit, the ion
temperature equation Eq. (A) reduces to

$\displaystyle\partial_{t}T_{i}+\frac{1}{\rho^{*}}\left[\phi,T_{i}\right]+U_{\parallel
i}\partial_{z}T_{i}$ $\displaystyle=+\frac{2}{3}T_{i}\left[\left(U_{\parallel
i}-U_{\parallel
e}\right)\frac{\partial_{z}n_{e}}{n_{e}}-\partial_{z}U_{\parallel e}\right]$
$\displaystyle+\partial_{z}\left(\chi_{\parallel
i}\partial_{z}T_{i}\right)+\frac{A_{\mathcal{E}}}{n_{e}}+\left(1-T_{i}\right)\frac{S_{\mathcal{N}}}{n_{e}},$
(33)

where the temperature anisotropy is neglected such that $T_{\parallel i}\simeq
T_{i}$ at high-collisionality Frei, Hoffmann, and Ricci (2022). Hence, Eq. (A)
is equivalent to the ion temperature equation given in Eq. (22b) used in the
Braginskii model where the source terms are due to the density and energy
sources (see Eq. (14)) considered in the ion GM hierarchy equation.

A high-collisional closure for the parallel ion heat flux can also be derived
using the Dougherty collision operator (or other collision operator models).
However, it has been shown that simplified collision operators can lead to
large deviations (with respect to the Fokker-Planck collision operator) in the
predictions of the conductivities Frei, Ernst, and Ricci (2022), especially
when using operators like Dougherty, as in this work.

Finally, it is worth noting that, while the Braginskii collisional closure
diverges in the low-collisionality limit Pitzal  _et al._ (2023), resulting in
an overestimate of the parallel heat flux, the GM hierarchy equation allows
for the evolution of $q_{\parallel i}$ using a larger number of GMs. However,
while this effect may not be as significant in LAPD due to the low
temperature, it could be more pronounced in the boundary region.

## References

  * Dorf and Dorr (2020) M. Dorf and M. Dorr, “Progress with the 5D full-F continuum gyrokinetic code COGENT,” Contributions to Plasma Physics 60, e201900113 (2020).
  * Hakim  _et al._ (2020) A. H. Hakim, N. R. Mandell, T. N. Bernard, M. Francisquez, G. W. Hammett, and E. L. Shi, “Continuum electromagnetic gyrokinetic simulations of turbulence in the tokamak scrape-off layer and laboratory devices,” Physics of Plasmas 27, 042304 (2020).
  * Michels  _et al._ (2022) D. Michels, P. Ulbl, W. Zholobenko, T. Body, A. Stegmeir, T. Eich, M. Griener, G. D. Conway, F. Jenko, and A. U. Team, “Full-f electromagnetic gyrokinetic turbulence simulations of the edge and scrape-off layer of ASDEX Upgrade with GENE-X,” Physics of Plasmas 29, 032307 (2022).
  * Ulbl  _et al._ (2023) P. Ulbl, T. Body, W. Zholobenko, A. Stegmeir, J. Pfennig, and F. Jenko, “Influence of collisions on the validation of global gyrokinetic simulations in the edge and scrape-off layer of TCV,” Physics of Plasmas 30 (2023).
  * Chang  _et al._ (2017) C. S. Chang, S. Ku, A. Loarte, V. Parail, F. Koechl, M. Romanelli, R. Maingi, J.-W. Ahn, T. Gray, J. Hughes, _et al._ , “Gyrokinetic projection of the divertor heat-flux width from present tokamaks to ITER,” Nuclear Fusion 57, 116023 (2017).
  * Stegmeir  _et al._ (2018) A. Stegmeir, D. Coster, A. Ross, O. Maj, K. Lackner, and E. Poli, “GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach,” Plasma Physics and Controlled Fusion 60, 035005 (2018).
  * De Oliveira  _et al._ (2022) D. S. De Oliveira, T. A. Body, D. Galassi, C. Theiler, E. Laribi, P. Tamain, A. Stegmeir, M. Giacomin, W. Zholobenko, P. Ricci, _et al._ , “Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case,” Nuclear Fusion 62, 096001 (2022).
  * Giacomin  _et al._ (2022) M. Giacomin, P. Ricci, A. Coroado, G. Fourestey, D. Galassi, E. Lanti, D. Mancini, N. Richart, L. N. Stenger, and N. Varini, “The gbs code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary,” Journal of Computational Physics , 111294 (2022).
  * Zeiler, Drake, and Rogers (1997) A. Zeiler, J. F. Drake, and B. Rogers, “Nonlinear reduced braginskii equations with ion thermal dynamics in toroidal plasma,” Physics of Plasmas 4, 2134 (1997).
  * Frei, Jorge, and Ricci (2020) B. J. Frei, R. Jorge, and P. Ricci, “A gyrokinetic model for the plasma periphery of tokamak devices,” Journal of Plasma Physics 86, 905860205 (2020).
  * Jorge  _et al._ (2019) R. Jorge, P. Ricci, S. Brunner, S. Gamba, V. Konovets, N. F. Loureiro, L. M. Perrone, and N. Teixeira, “Linear theory of electron-plasma waves at arbitrary collisionality,” Journal of Plasma Physics 85, 905850211 (2019).
  * Frei, Hoffmann, and Ricci (2022) B. J. Frei, A. C. D. Hoffmann, and P. Ricci, “Local gyrokinetic collisional theory of the ion-temperature gradient mode,” Journal of Plasma Physics 88, 905880304 (2022).
  * Frei  _et al._ (2023) B. J. Frei, A. C. D. Hoffmann, P. Ricci, S. Brunner, and Z. Tecchioll, “Moment-based approach to the flux-tube linear gyrokinetic model,” Journal of Plasma Physics 89, 905890414 (2023).
  * Jorge, Frei, and Ricci (2019) R. Jorge, B. J. Frei, and P. Ricci, “Nonlinear gyrokinetic Coulomb collision operator,” Journal of Plasma Physics 85, 905850604 (2019).
  * Frei  _et al._ (2021) B. J. Frei, J. Ball, A. C. D. Hoffmann, R. Jorge, P. Ricci, and L. Stenger, “Development of advanced linearized gyrokinetic collision operators using a moment approach,” Journal of Plasma Physics 87, 905870501 (2021).
  * Frei, Ernst, and Ricci (2022) B. J. Frei, S. Ernst, and P. Ricci, “Numerical implementation of the improved sugama collision operator using a moment approach,” Physics of Plasmas 29, 093902 (2022).
  * Hoffmann, Frei, and Ricci (2023a) A. C. D. Hoffmann, B. J. Frei, and P. Ricci, “Gyrokinetic simulations of plasma turbulence in a z-pinch using a moment-based approach and advanced collision operators,” Journal of Plasma Physics 89, 905890214 (2023a).
  * Hammett and Perkins (1990) G. W. Hammett and F. W. Perkins, “Fluid moment models for Landau damping with application to the ion-temperature-gradient instability,” Physical Review Letters 64, 3019 (1990).
  * Beer and Hammett (1996) M. A. Beer and G. W. Hammett, “Toroidal gyrofluid equations for simulations of tokamak turbulence,” Physics of Plasmas 3, 4046 (1996).
  * Brizard and Hahm (2007) A. J. Brizard and T. S. Hahm, “Foundations of nonlinear gyrokinetic theory,” Reviews of Modern Physics 79, 421 (2007).
  * Hoffmann, Frei, and Ricci (2023b) A. C. D. Hoffmann, B. J. Frei, and P. Ricci, “Dimits shift and impact of collisions using a gyrokinetic moment-based approach,” arXiv preprint arXiv:2308.01016 (2023b).
  * Mandell, Dorland, and Landreman (2018) N. R. Mandell, W. Dorland, and M. Landreman, “Laguerre–Hermite pseudo-spectral velocity formulation of gyrokinetics,” Journal of Plasma Physics 84, 905840108 (2018).
  * Mandell  _et al._ (2022) N. Mandell, W. Dorland, I. Abel, R. Gaur, P. Kim, M. Martin, and T. Qian, “GX: a GPU-native gyrokinetic turbulence code for tokamaks and stellarators,” arXiv:2209.06731 (2022).
  * Gekelman  _et al._ (1991) W. Gekelman, H. Pfister, Z. Lucky, J. Bamber, D. Leneman, and J. Maggs, “Design, construction, and properties of the large plasma research device- The LAPD at UCLA,” Review of scientific instruments 62, 2875 (1991).
  * Gilmore  _et al._ (2015) M. Gilmore, A. G. Lynn, T. R. Desjardins, Y. Zhang, C. Watts, S. C. Hsu, S. Betts, R. Kelly, and E. Schamiloglu, “The helcat basic plasma science device,” Journal of Plasma Physics 81, 345810104 (2015).
  * Furno  _et al._ (2017) I. Furno, R. Agnello, U. Fantz, A. Howling, R. Jacquier, C. Marini, G. Plyushchev, P. Guittienne, and A. Simonin, “Helicon wave-generated plasmas for negative ion beams for fusion,” in  _EPJ Web of Conferences_ , Vol. 157 (EDP Sciences, 2017) p. 03014.
  * Tynan  _et al._ (2004) G. R. Tynan, M. J. Burin, C. Holland, G. Antar, N. Crocker, and P. H. Diamond, “Radially sheared azimuthal flows and turbulent transport in a cylindrical plasma,” Physics of plasmas 11, 5195 (2004).
  * Castellanos  _et al._ (2005) O. F. Castellanos, E. Anabitarte, J. M. Senties, C. Hidalgo, and M. A. Pedrosa, “Parallel flows and turbulence in a linear plasma machine,” Plasma Physics and Controlled Fusion 47, 2067 (2005).
  * Rogers and Ricci (2010) B. N. Rogers and P. Ricci, “Low-frequency turbulence in a linear magnetized plasma,” Physical Review Letters 104, 225002 (2010).
  * Chodura (1982) R. Chodura, “Plasma–wall transition in an oblique magnetic field,” The Physics of Fluids 25, 1628 (1982).
  * Loizu  _et al._ (2012) J. Loizu, P. Ricci, F. D. Halpern, and S. Jolliet, “Boundary conditions for plasma fluid models at the magnetic presheath entrance,” Physics of Plasmas 19, 122307 (2012).
  * Geraldini, Parra, and Militello (2017) A. Geraldini, F. I. Parra, and F. Militello, “Gyrokinetic treatment of a grazing angle magnetic presheath,” Plasma Physics and Controlled Fusion 59, 025015 (2017).
  * Friedman  _et al._ (2012) B. Friedman, T. A. Carter, M. V. Umansky, D. Schaffner, and B. Dudson, “Energy dynamics in a simulation of lapd turbulence,” Physics of Plasmas 19, 102307 (2012).
  * Popovich  _et al._ (2010) P. Popovich, M. V. Umansky, T. A. Carter, and B. Friedman, “Modeling of plasma turbulence and transport in the large plasma device,” Physics of Plasmas 17, 122312 (2010).
  * Fisher  _et al._ (2015) D. M. Fisher, B. N. Rogers, G. D. Rossi, D. S. Guice, and T. A. Carter, “Three-dimensional two-fluid braginskii simulations of the large plasma device,” Physics of Plasmas 22, 092121 (2015).
  * Ross  _et al._ (2019) A. Ross, A. Stegmeir, P. Manz, D. Groselj, W. Zholobenko, D. Coster, and F. Jenko, “On the nature of blob propagation and generation in the large plasma device: Global grillix studies,” Physics of Plasmas 26, 102308 (2019).
  * Shi  _et al._ (2017) E. L. Shi, G. W. Hammett, T. Stoltzfus-Dueck, and A. Hakim, “Gyrokinetic continuum simulation of turbulence in a straight open-field-line plasma,” Journal of Plasma Physics 83, 905830304 (2017).
  * Pan  _et al._ (2018) Q. Pan, D. Told, E. Shi, G. W. Hammett, and F. Jenko, “Full-f version of GENE for turbulence in open-field-line systems,” Physics of Plasmas 062303, 1 (2018).
  * Dougherty (1964) J. P. Dougherty, “Model Fokker-Planck equation for a plasma and its solution,” The Physics of Fluids 7, 1788 (1964).
  * Carter (2006) T. A. Carter, “Intermittent turbulence and turbulent structures in a linear magnetized plasma,” Physics of Plasmas 13, 010701 (2006).
  * Brizard (1992) A. J. Brizard, “Nonlinear gyrofluid description of turbulent magnetized plasmas,” Physics of Fluids B: Plasma Physics 4, 1213 (1992).
  * Snyder and Hammett (2001) P. B. Snyder and G. W. Hammett, “A landau fluid model for electromagnetic plasma microturbulence,” Physics of Plasmas 8, 3199 (2001).
  * Strintzi and Scott (2004) D. Strintzi and B. Scott, “Nonlocal nonlinear electrostatic gyrofluid equations,” Physics of plasmas 11, 5452 (2004).
  * Madsen (2013) J. Madsen, “Full-f gyrofluid model,” Physics of Plasmas 20, 072301 (2013).
  * Held  _et al._ (2016) M. Held, M. Wiesenberger, J. Madsen, and A. Kendl, “The influence of temperature dynamics and dynamic finite ion larmor radius effects on seeded high amplitude plasma blobs,” Nuclear Fusion 56, 126005 (2016).
  * Wiesenberger  _et al._ (2019) M. Wiesenberger, L. Einkemmer, M. Held, A. Gutierrez-Milla, X. Sáez, and R. Iakymchuk, “Reproducibility, accuracy and performance of the feltor code and library on parallel computer architectures,” Computer Physics Communications 238, 145 (2019).
  * Galassi  _et al._ (2022) D. Galassi, C. Theiler, T. Body, F. Manke, P. Micheletti, J. Omotani, M. Wiesenberger, M. Baquero-Ruiz, I. Furno, M. Giacomin, _et al._ , “Validation of edge turbulence codes in a magnetic x-point scenario in torpex,” Physics of Plasmas 29 (2022).
  * Wiesenberger and Held (2022) M. Wiesenberger and M. Held, “Long-wavelength closures for collisional and neutral interaction terms in gyro-fluid models,” Journal of Physics: Conference Series 2397, 012015 (2022).
  * Wiesenberger and Held (2023) M. Wiesenberger and M. Held, “Effects of plasma resistivity in three-dimensional full-f gyro-fluid turbulence simulations,” arXiv preprint arXiv:2309.10414 (2023).
  * Held, Wiesenberger, and Kendl (2020) M. Held, M. Wiesenberger, and A. Kendl, “Padé-based arbitrary wavelength polarization closures for full-F gyro-kinetic and-fluid models,” Nuclear Fusion 60, 066014 (2020).
  * Held and Wiesenberger (2022) M. Held and M. Wiesenberger, “Beyond the oberbeck–boussinesq and long wavelength approximation,” Nuclear Fusion 63, 026008 (2022).
  * Jorge, Ricci, and Loureiro (2017) R. Jorge, P. Ricci, and N. F. Loureiro, “A drift-kinetic analytical model for scrape-off layer plasma dynamics at arbitrary collisionality,” Journal of Plasma Physics 83, 905830606 (2017).
  * Goerler  _et al._ (2011) T. Goerler, X. Lapillonne, S. Brunner, T. Dannert, F. Jenko, F. Merz, and D. Told, “The global version of the gyrokinetic turbulence code GENE,” Journal of Computational Physics 230, 7053 (2011).
  * Sarazin  _et al._ (2010) Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, P. Ghendrih, G. Latu, A. Strugarek, and G. Dif-Pradalier, “Large scale dynamics in flux driven gyrokinetic turbulence,” Nuclear Fusion 50, 054004 (2010).
  * Tripathi, Pribyl, and Gekelman (2011) S. K. P. Tripathi, P. Pribyl, and W. Gekelman, “Development of a radio-frequency ion beam source for fast-ion studies on the large plasma device,” Review of Scientific Instruments 82, 093501 (2011).
  * Gradshteyn and Ryzhik (2014) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).
  * Wiesenberger and Held (2020) M. Wiesenberger and M. Held, “Angular momentum and rotational energy of mean flows in toroidal magnetic fields,” Nuclear Fusion 60, 096018 (2020).
  * Yu, Krasheninnikov, and Guzdar (2006) G. Q. Yu, S. I. Krasheninnikov, and P. N. Guzdar, “Two-dimensional modelling of blob dynamics in tokamak edge plasmas,” Physics of Plasmas 13, 042508 (2006).
  * Bodi  _et al._ (2011) K. Bodi, G. Ciraolo, F. Schwander, and E. Serre, “Impact of the boussinesq approximation in tokamak scrape-off layer turbulence,” in  _38th Conference on Plasma Physics_ , Vol. 1 (2011) p. 121.
  * Paruta  _et al._ (2018) P. Paruta, P. Ricci, F. Riva, C. Wersal, C. Beadle, and B. Frei, “Simulation of plasma turbulence in the periphery of diverted tokamak by using the GBS code,” Physics of Plasmas 25, 112301 (2018).
  * Stangeby (1995) P. C. Stangeby, “The bohm–chodura plasma sheath criterion,” Physics of Plasmas 2, 702 (1995).
  * Mosetto  _et al._ (2015) A. Mosetto, F. D. Halpern, S. Jolliet, J. Loizu, and P. Ricci, “Finite ion temperature effects on scrape-off layer turbulence,” Physics of Plasmas 22, 012308 (2015).
  * Riemann (1981) K.-U. Riemann, “Kinetic theory of the plasma sheath transition in a weakly ionized plasma,” The Physics of Fluids 24, 2163 (1981).
  * Baalrud  _et al._ (2015) S. D. Baalrud, B. Scheiner, B. Yee, M. Hopkins, and E. Barnat, “Extensions and applications of the Bohm criterion,” Plasma Physics and Controlled Fusion 57, 044003 (2015).
  * Churchill  _et al._ (2017) R. M. Churchill, J. M. Canik, C. S. Chang, R. Hager, A. W. Leonard, R. Maingi, R. Nazikian, and D. P. Stotler, “Kinetic simulations of scrape-off layer physics in the DIII-D tokamak,” Nuclear Materials and Energy 12, 978 (2017).
  * Arakawa (1997) A. Arakawa, “Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I,” Journal of computational physics 135, 103 (1997).
  * Parker  _et al._ (1993) S. E. Parker, R. J. Procassini, C. K. Birdsall, and B. I. Cohen, “A suitable boundary condition for bounded plasma simulation without sheath resolution,” Journal of Computational Physics 104, 41 (1993).
  * Gekelman  _et al._ (2016) W. Gekelman, P. Pribyl, Z. Lucky, M. Drandell, D. Leneman, J. Maggs, S. Vincena, B. Van Compernolle, S. K. P. Tripathi, G. Morales, _et al._ , “The upgraded large plasma device, a machine for studying frontier basic plasma physics,” Review of Scientific Instruments 87, 025105 (2016).
  * Chapman and Cowling (1941) S. Chapman and T. G. Cowling, “The velocity of diffusion in a mixed gas; the second approximation,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 179, 159 (1941).
  * Braginskii (1965) S. I. Braginskii, “Transport processes in a plasma,” Reviews of plasma physics 1, 205 (1965).
  * Pitzal  _et al._ (2023) C. Pitzal, A. Stegmeir, W. Zholobenko, K. Zhang, and F. Jenko, “Landau-fluid simulations of edge-SOL turbulence with GRILLIX,” Physics of Plasmas 30, 122502 (2023).