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

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.


I. INTRODUCTION
Despite recent progress in the development of gyrokinetic (GK) codes, such as COGENT 1 , Gkeyll 2 , GENE-X 3,4 and XGC 5 , 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.6-8) based on the drift-reduced Braginskii model 9 , 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. 10 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 [11][12][13] , which are absent in Braginskii-like fluid models, and collisional effects modeled using advanced collision operators [13][14][15][16][17] 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 18,19 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 δ f regime (where only the a priori small deviation of the distribution function from thermal equilibrium is evolved 20 ) showing convergence with a low number of GMs, in particular at high collisionality 13 .We notice that the first nonlinear δ f simulations using advanced linearized collision operators have been recently performed 17,21 and that a similar approach has also been implemented to perform δ f turbulent calculations focusing on the core region 22,23 .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 24 , HelCat 25 and RAID 26 , 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 24,[26][27][28][29] .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 [30][31][32] and by the low plasma temperatures comparable to typical SOL values (e.g., T i ≲ T e ∼ 6 eV in typical LAPD discharges 33 ), 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 29,[34][35][36] , 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 velocityspace 37 .LAPD turbulent simulations using the GK GENE code are also reported in Ref. 38 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. 10. Consistently with Refs.37  and 38, 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 39 collision operator (similar to the one used in Refs.37 and 38).In addition, we neglect FLR effects and we refer to our ion model as a GK model for consistency with previous works 37,38 , 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 9 .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 37,38 , 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 31 are implemented in the parallel direction, which are equivalent to the ones used in the previous Braginskii simulation of LAPD 29,35 .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 9 (or simply Braginskii model), similarly to Refs.29 and 35, 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. 29, 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 35,40 and with previous GK simulations 37,38 falls outside the scope of our study, but will be addressed in future work.
We remark that the GM approach of Ref. 10 shares some similarities with previous (full-F or δ f ) gyrofluid models 19,[41][42][43][44][45][46][47][48][49] , 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 [49][50][51] and collisional closures 48,52 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 ∥e , and temperature T e .A vorticity equation is derived using the Boussinesq approximation for the electrostatic potential φ , 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 crosssection, we define the cartesian coordinate system (x, y, z), such that the (x, y) coordinates describe the plane perpendicular to 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 B = ∇ × A = Bz, where A is the magnetic vector potential, which is constant in time, and 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 A and we derive the ion full-F GM hierarchy equation in Sec.II B. A presentation of the reduced cold-ion model is then obtained from the GM hierarchy equation in Sec.II C. The fluid electron model follows in Sec.II D and the vorticity equation is derived in Sec.II E. Sec.II F describes the two-fluid Braginskii model used for comparison purposes and, finally, Sec.II G details the Bohm sheath boundary conditions we use in our simulations.

A. Ion full-F model
Focusing on the electrostatic and neglecting FLR effects with constant and straight magnetic field lines, the ion oneform Γ i 10 expressed in the gyrocenter coordinates Z = (R g , µ, v ∥ , θ ), where R g is the gyrocenter position, µ = m i v 2 ⊥ /(2B) is the magnetic moment and v ∥ = b • v is the velocity parallel to the magnetic field with b = B/B, reduces to with q i Φ i = q i φ + µB and q i A * i = q i A + m i v ∥ b.We remark that, in Eq. ( 1), the electrostatic potential φ is evaluated at the gyrocenter position, i.e. φ = φ (R g ), such that ion FLR effects are neglected.From Eq. ( 1), we deduce the ion equations of motion, and μ = 0, with v ∥ = Ṙg • b and E = −∇φ the electric field.Eq. (2a) describes the parallel streaming along the magnetic field lines and the perpendicular drift due to the E × B velocity, while Eq.(2b) represents the acceleration in the parallel direction associated with the electric field E. Using the equations of motion given in Eq. ( 2), the evolution equation of the full-F (gyrophase-independent) ion distribution function, ), in the electrostatic limit of the GK ion Boltzmann equation 20 is given by where J i = B/m i is the gyrocenter phase-space Jacobian, which is a constant in the case of linear devices.On the righthand side of Eq. ( 3), S i = S i (R g , µ, v ∥ ) = S N + S E models the particle (S N ) and energy (S E ) sources and are defined by 53,54 respectively.We note that the sources, given in Eq. ( 4), neglect FLR and polarization effects associated with the transformation from particle to gyrocenter coordinate 48 .In Eq. ( 4), the perpendicular normalized velocity-space coordinate is defined by x i = µB/T i0 , while we use 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 ≪ T e ) is used and the ion energy source is neglected.
The functions A N = A N (x, y) and A E = A 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 55 .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 tophat-like shape.For instance, A N (x, y) is given by 29 where r = 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, A N0 is a positive and constant coefficient, which represents the physical particle fuelling rate near the center of the device, while A N ∞ represents a small (A N ∞ ≪ A N 0 ) positive and constant particle source away from r ∼ r s added for numerical reasons, in particular, to avoid regions of small or negative plasma density.Similar definitions for A E (x, y), A E 0 and A E ∞ are used.In Eq. ( 4), we also introduce a shifted Maxwellian distribution function defined by 41,44 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 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 39 , given by where T i = dvF i m i (v − u i ) 2 /(3N i ) and u i = dvF i v/N i are the ion temperature and mean fluid velocity, and ) is the ion-ion collision frequency, which is constant in the present work.The effects of ionelectron collisions are neglected in Eq. (3) since they occur on a time scale larger by, at least, a factor proportional to 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 37 and in the GENE 38 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 ∥ , µ), 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.

B. Full-F ion GM Hierarchy Equation
Following Ref. 10, 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 ∥i ), and Laguerre, L j (x i ), velocity-space polynomials 56 , such that where N p j are the ion GMs, evaluated by using the Hermite and Laguerre orthogonality relations 56 By introducing the GM projector 10,22,44,52 with χ = χ(R g , µ, v ∥ ,t) being an arbitrary gyrocenter phasespace function, we find N p j = ∥1∥ p j from Eq. ( 10).We remark that, in Eq. ( 8), the shift of the parallel velocity coordinate s ∥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 C).These flows appear at the sheath entrance where ions are accelerated to the ion sound speed (see Sec. II G).Additionally, we note that F Mi , defined in Eq. ( 6), is assumed to have the same parallel and perpendicular temperature, T ∥i = T ⊥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 ≲ 1 eV) compared to the boundary region in fusion devices (where T i ≳ 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 ∥i , and the ion parallel and perpendicular pressure and temperature P ∥i = T ∥i N i and P ⊥i = T ⊥i N i , respectively.Indeed, using Eq. ( 9), we derive that with the total ion temperature defined by We remark that Eq. (11b) is a direct consequence of our choice of using a shifted parallel velocity-space coordinate s ∥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 14,48,52 , 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, N p j .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 = 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 B), the potential φ to T e0 /e, the parallel and perpendicular spatial scales to R and ρ s0 = c s0 /Ω 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 ∥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 N p j , i.e.
where the GM projections are given by p with ρ * = ρ s0 /R and τ i = T i0 /T e0 .In Eq. ( 13), we introduce the Poisson bracket operator that is The GM expansions of the particle and energy sources, S p j N and S p j E , are given by 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 ∥i , x i ) and project it onto the Hermite-Laguerre basis.This yields where T i is expressed in terms of the GMs using Eq. ( 11).The nonlinear Dougherty collision operator conserves particles (C 00 i = 0), momentum (C 10 i = 0) and energy (C 20 i = √ 2C 01 i ).While simpler in form compared to the GM expansion of the nonlinear Fokker-Planck Landau collision operator 14 , 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 14 will be considered in future work.
To obtain the time evolution of the GMs N p j , it is necessary to derive an explicit expression for the time derivative of the ion parallel velocity, ∂ t U ∥i which appears in Eq. ( 12) and resulting from the use of the shifted parallel velocity-space coordinate s ∥i .By setting (p, j) = (1, 0) in Eq. ( 12) and using the fact that N 10 vanishes exactly (see Eq. (11b)), we derive the desired expression for ∂ t U ∥i given by where the parallel ion pressure P ∥i is expressed in terms of GMs according to Eq. (11c).We note that the full-F GM hierarchy equation, given in Eq. ( 12), can also be derived from the electromagnetic full-F GM hierarchy equation described in Ref. 10.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 ∥i and T ⊥i in Eq. (11c) and Eq.(11d), respectively) as in Ref. 10, additional terms proportional to their time derivatives and gradients appear in Eq. (12).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. ( 12) 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. ( 15)) as our model neglects FLR effects yielding additional coupling in j 12 .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 N p j = 0 for all (p, j) > (P, J) with 0 ≤ P, J < ∞.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.

C. 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.29, 34, and 35) where the effects of finite ion temperature T i are neglected.
In the cold ion limit (T i = 0), only the GM N 00 i and U ∥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 ∥i in Eq. ( 16) can be neglected.As a consequence, the ion GM hierarchy equation given in Eq. ( 12) reduces to the ion gyrocenter continuity equation for N i and to the ion parallel momentum equation for U ∥i 47,48,57 , i.e.
respectively.We remark that the particle and momentum conservation of the collision operator is used in deriving Eq. ( 17).

D. Electron fluid model
We use the Braginskii model to evolve the electron dynamics, avoiding the evolution of their distribution function, in contrast to Refs.37 and 38.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 ⊥ ρ s ∼ 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 ∥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 where the normalized parallel electrical resistivity and electron thermal conductivity are given by ν ∥ = ν 0 /T 3/2 e and χ ∥e = 1.075T 5/2 e /ν 0 , respectively.Here, .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 ∥ = n e (U ∥i −U ∥e ).

E. Vorticity equation
We now obtain the vorticity equation that governs the evolution of the electrostatic potential φ .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 10,20 −en e + q i N i = −∇ • 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 51 , requires the development of FLR closures at arbitrary wavelength (see, e.g., Ref. 50) 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 37,38 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 8,36 and is also considered below to derive the vorticity equation.While Eq. ( 19) can be solved to obtain φ 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 29,[34][35][36] .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 where we introduce the vorticity variable Ω = ∇ 2 ⊥ φ .The vorticity equation is then We note that the effects of the Boussinesq approximation on plasma turbulence are the subject of previous studies 36,51,[58][59][60] .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 φ from the vorticity variable Ω.We use the vorticity equation given in Eq. ( 21) to evolve Ω when considering the full-F ion GM hierarchy equation and the cold ion models, given in Eqs. ( 12) and ( 17) respectively, coupled to the fluid electron model in Eq. ( 18).

F. Two-fluid Braginskii fluid model
We finally introduce the two-fluid Braginskii fluid model 9 , 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 ≃ 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 ∥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 ∥e and T e already described in Sec.II D, the two-fluid Braginskii equations prescribe a parallel ion momentum equation to evolve U ∥i , an ion temperature equation to evolve T i , and vorticity equations for Ω.These equations are given by respectively.In Eq. (22b), χ ∥i = 1.32 m e /m i (τ i T i ) 5/2 /ν 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 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 ∂ z φ , appearing in the parallel ion momentum equation, Eq. (17b), is approximated in Eq. (22a) by the electron parallel pressure gradient, such that ∂ z φ ≃ ∂ 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. τ i ∇ 2 ⊥ T i , present in the vorticity equation, Eq. (22c), 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.

G. Boundary conditions
Boundary conditions are required for the ion GMs, N p j , the electron fluid quantities, N e , U ∥e , T e , and the potential φ in the perpendicular (x, y) plane at x = ±L x /2 and y = ±L y /2 and at the end plates located in the z direction at z = ±L z /2, where a sheath forms due to the plasma-wall interaction.
At x = ±L x /2 and y = ±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 61 .Since the sheath region cannot be modeled by the field equations derived in Sec.II E (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 37,38 , 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 61 which sets the value of the parallel electron and ion velocities, U ∥e and U ∥i , at the sheath entrance.Therefore, we assume that 31,62 U ∥e (x, y, z = ±L z /2) = ± T e,s e Λ−φ s /T e,s , (23a) with Λ = log m i /(2m e ) ≃ 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 = ±L z /2) and T i,s = T i (x, y, z = ±L z /2), and, similarly φ s = φ (x, y, z = ±L z /2).
We notice that the boundary conditions in Eq. ( 23) reduce to the ones used in Ref. 29 when T i ≪ T e and correspond to the ones used in SOL turbulent simulations using the driftreduced Braginskii model 8 .For the remaining quantities, we assume, for simplicity, that the gradients of electron density, n e , electron temperature, T e , ion GMs, N p j , and electrostatic potentials, φ , vanish along the direction of the magnetic field at the sheath entrance, i.e. homogenous Neumann boundary conditions are imposed at z = ±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.31 and 62, can be extended to an arbitrary number of GMs and kinetic sheath boundary conditions can also be developed [63][64][65] .Magnetic field lines intercept the machine wall with a small oblique angle in fusion devices, further complicating the treatment of the sheath boundary conditions 30,32 .

III. NUMERICAL IMPLEMENTATION
To solve the full-F ion GM hierarchy given in Eq. ( 12) coupled with the electron fluid model in Eq. ( 18) and the vorticity equation in Eq. ( 21), we have developed a new threedimensional 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 C) and the two-fluid Braginskii model (Sec.II F) for comparison with the GM results.
To evolve the plasma dynamics, we employ similar numerical algorithms as the two-fluid GBS code 8 .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 is evaluated by using a fourth-order Arakawa method 66 .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 8 .To avoid checkerboard patterns 60 , the grid used to evolve the parallel velocities, U ∥e and U ∥i , and the GMs N p j with odd p, is staggered to the left along the z-direction by ∆z/2 (∆z is the grid spacing) with respect to the grid where the other fluid quantities, i.e. n e , T e , Ω (and thus φ ), and the GMs N p j with even p are evaluated.Fourth-order interpolation techniques are used between two staggered grids 60 .To improve the numerical stability of our numerical simulations, parallel and perpendicular numerical diffusions, such as 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, η ⊥ and η 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 = 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 = A 0 [1 − tanh {(r − r s )/(4L s )}]/2.We use A 0 = 0.4 for f = n e and A 0 = 0.4 for f = T e .In addition, we set φ = ΛT e to avoid unphysical and large electron current into the sheath region.The initial values of N 20 and N 01 , given the initial ion density N and ion temperature T i profiles, are obtained by inverting Eq. (11), which yields Finally, the parallel velocities, U ∥i and U ∥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 N 01 , N 02 ), and parallel velocities, U ∥e and U ∥i , profiles to seed turbulence.Typically, a quasi-steady state is achieved after 100 c s0 /R time unit (corresponding to t ∼ 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 B), with different numbers of GMs and values of collisionality and compare them with the Braginskii model introduced in Sec.II F.
Our simulations parameters are similar to those used in Ref. 29, where a helium LAPD plasma is considered.These parameters are sumarized as follows: n e0 = 2 × 10 12 cm −3 , T e0 = 6 eV, T i0 = 3 eV (τ i = 0.5), Ω i ∼ 960 kHz, ρ s0 = 1.4 cm, c s0 = 1.3 × 10 6 cm s −1 , m i /m e = 400 and ν 0 = 0.03.The LAPD vacuum chamber has a radius R ≃ 0.56 m (i.e., R ≃ 40ρ s0 ) and a parallel length of L z ≃ 18 m, such that we use L x = L y = 100ρ s0 (or L x ∼ L y ∼ 1.4 m) and L z = 36R.The reference time is R/c s0 ∼ 43 µs.We consider the following parameters for the density and temperature sources L s = 1ρ s0 , r s = 20ρ s0 , A N 0 = A T e0 = 0.04 (with A N ∞ = A T e ∞ = 0.001) and A E 0 = 0.02 (with A E ∞ = A N ∞ ).We use a numerical resolution of N x = N y = 192 in the perpendicular plane, which corresponds to a grid spacing of ∼ 0.52ρ s0 .This resolution is also sufficient to resolve turbulence occurring at perpendicular wavelengths longer than ρ s0 in the perpendicular plane, which is dominated by the long-wavelength KH instability 29 .On the other hand, a coarser resolution of N z = 64 is used in the parallel direction.Given our numerical resolution, we use η ⊥ = 0.05 and η 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 ν 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 34.In this regime, the ion mean-free-path, λ mp f , is considerably shorter than the total length L z , i.e. λ mp f /L z ≃ √ 2τ i R/L z /ν i ≪ 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 ν i ≃ 4 × 10 −3 yielding λ mp f /L z ∼ 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. ( 15) 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 12 and trapped particles 13 (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 10 .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 ≥ 2 GMs is from the pitch-angle scattering contained in the collision operator given in Eq. ( 15), which is dominated by collisional damping.Indeed, if collisions are neglected GMs with j > 1 are completely decoupled from the the GMs with j ≤ 1 This section is structured as follows.First, Sec.IV A 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 B. Finally, we investigate the ion distribution function in velocity-space in Sec.IV C and the GM spectrum in quasi-steady state in Sec.IV D as a function of the number of GMs and for the two collisionality regimes.

A. 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. ( 12), 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 ≲ 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 ∼ 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 ⊥ ρ s0 ∼ 0.5 (k ⊥ 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 34 , nonlinearly trigger a KH instability, characterized by a long perpendicular wavelength (k ⊥ ρ s0 ≪ 1) and k ∥ ≃ 0. The KH instability becomes clearly visible around t ≃ 45R/c s0 .This instability, which has been shown to dominate the radial transport in LAPD 29 , saturates at t ≃ 53R/c s0 , transporting the plasma to the r ≳ 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. 29) when the nonlinear term [φ , Ω] in Eq. ( 21) is artificially suppressed.
After t ∼ 130R/c s0 (similarly to previous GK and Braginskii turbulent simulations of LAPD 34,38 ), 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 ∼ 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 timeaveraged 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 ∼ 41500 CPU hours on the multi-core partition of the Piz Daint supercomputer, approximatively.
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 φ , electron density n e , and electron temperature T e reveal elongated structures.All quantities show weak parallel gradients k ∥ ≃ 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 29,34,35 and GK 37,38 turbulent simulations of LAPD.
We now examine the time-averaged radial profiles.These profiles are obtained by averaging over a time window of ∼ 2 ms during the quasi-steady state as well as over the central region of LAPD, i.e. −8R ≤ z ≤ +8R (or −4m ≲ z ≲ 4m), a region commonly considered to present experimental data 40 (a similar approach is used in previous GK simulations 37,38 ).More precisely, the averaged profile of the electrostatic potential φ is computed by defining ⟨φ ⟩ t,z = +ℓ z /2 −ℓ z /2 dz τ 0 dtφ /(ℓ z τ), with ℓ z and τ are the z and time length of the average windows, respectively.The results are shown in Fig. 4, which displays the profiles of ⟨φ ⟩ t,z , ⟨n e ⟩ t,z , and ⟨T e ⟩ 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 29 .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 ∼ r s , where the fluctuation level is large (see Sec. IV B).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 dif-ferent 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 29 and GK simulations 37,38 .
We remark that the electrostatic potential profile φ follows approximatively the electron temperature T e , as shown in Fig. 4. Indeed, φ ∼ ΛT e is required to have comparable electron and ion outflows in steady-state, such that U ∥i ∼ U ∥e near the end plates, according to Eq. ( 23).To verify that φ ∼ ΛT e in our simulations, we evaluate the radial profile of the instantaneous difference, φ − Λ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 quasisteady state.The results are shown in Fig. 5.We first observe that the GM and Braginskii simulations yield similar negative φ − ΛT e values.On the other hand, φ − ΛT e is roughly constant and approximatively vanishes for all radii in the coldion model.The deviations observed in φ − Λ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 φ −ΛT e ≃ −T e,s ln (1 + τ i T i,s /T e,s ) /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 φ − ΛT e ≃ 0 in the cold-ion model.Consequently, due to the finite ion temperature dependence of the boundary condition Eq. ( 23b), the equilibrium electrostatic potential φ exhibits higher amplitude in the cold ion model, which is confirmed in the timeaveraged profiles observed in Eq. ( 4).Finally, in all cases considered, φ − ΛT e remains smaller than the values of φ and ΛT e (φ − ΛT e ∼ 0.1 for r ≲ r s compared to φ ∼ ΛT e ∼ 2, see Fig. 3).

B. Turbulence analysis
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, φ , is defined by φ = φ − ⟨φ ⟩ t , where ⟨φ ⟩ t = τ 0 dtφ /τ denotes the time-averaged potential.Similar definitions for the other quantities are used.The top panels of Fig. 6 show instantaneous snapshots of φ in the plane perpendicular at the center of the device z = 0R, while the bottom panels illustrate φ snapshots.The Braginskii, cold-ion, and GM simulations with various (P, J) are considered.
We first observe that the fluctuations in the Braginskii model exhibits similar structures than those obtained in Ref. 35.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 ∼ r s (see Fig. 9).Notably, the φ 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 φ 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 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 φ 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 φ /T e ≲ 1 at all radii with, for example, φ ∼ 0.5 and T e ∼ 0.6 for r ≲ r s .
In Fig. 8, we present snapshots of the vorticity for the E × B flow, which is described by the vorticity variable Ω = ∇ 2 ⊥ φ , 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 ⊥ ρ s0 ≪ 1, driven by the KH instability, which develop in the region of steep equilibrium gradients near r ∼ 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.
We now proceed to analyze the root mean square (RMS) of the fluctuations, defined as n e 2 t in the case of the electron density fluctuations n e (and similarly for the other quantities).Fig. 9 displays the RMS of the electron density, n e , and the electrostatic potential, φ , fluctuations plotted as a function of the radius.The data are computed at z = 0R and normalized to ⟨n e ⟩ t (r) (and to ⟨φ ⟩ t (r)) 33,37 .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 35,36 and GK simulations 37 .Consistent with the observations made in Fig. 7, the RMS values reach their maximum when the gradients are most pronounced near r ∼ r s .For r ≲ r s and for r ≳ 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 φ ).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 n e ) are obtained in the cold-ion model.We compare the RMS of the parallel electrical current J ∥ 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 J ∥ 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 ∥ = 0 is imposed everywhere 67 .We remark that larger fluctuations of J ∥ 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.
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 n e .The skewness of the density is often used to characterize the dominance of plasma holes and blobs, associated with negative and positive skewness respectively 33,35,36 .Fig. 11 shows the skewness of the electron density n e .In all cases, the skewness is negative for r ≲ 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 ≳ 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 36 and GK 37,38 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 35,36 and validated with experimental data 33 .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.

C. Ion distribution function at quasi-steady state
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 ∥ /v Ti = s ∥i + √ 2τ i U ∥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 = ±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 ∥ ≃ 0 at z = 0R, where U ∥i ≃ 0. The absence of fine velocity-space structures in Fig. 12 in both the HC and LC regimes confirms the weak dependence of turbulence prop- erties on the number of GMs, as reported in Sec.IV B. 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. 38).The Maxwellian shape of F i observed in our simulations is also consistent with the fact that both two-fluid Braginskii 29,36 and GK 37,38 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 ∥ /v Ti ∼ 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. ( 12) does not enforce the positivity of the ion distribution function when evolving a finite number of GMs.For instance, the use of v ∥ /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 ∥ /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 N p j v ∥ , defined with respect to v ∥ /v Ti as the argument of H p , i.e.
are used to expand F i , it is found that N p j v ∥ ̸ = 0 for (p, j) > 0, even when F i is a Maxwellian distribution function centered at U ∥i ̸ = 0. Indeed, using Eq. ( 25), one derives the analytical expression of the unshifted GMs for where U ∥i is normalized to c s0 .While the amplitude of the unshifted GM decreases rapidly in the presence of subsonic ion flow, U ∥i ≪ 1, the decrease of the amplitude with p is slower in the presence of sonic flows, such that N p j v ∥ ∼ 2 p /p!.

D. GM spectrum at quasi-steady sate
To better assess the velocity-space representation of F i in our simulations, we plot the absolute value of the GMs, N p j (normalized to 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 12 , the amplitude of the GMs is expected to follow N p j ∼ exp [−c p (ν i )p α ] in the p ≫ 1 limit (with c p (ν i ) a positive coefficient that increases with collisionality, ν i and α > 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 velocityspace resolution of F i .On the other hand, the contributions from N p0 with p ≳ 4 are negligible in the HC regime, thereby justifying the closure by truncation for P ≳ 4. We also notice that 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 B 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 11,12 , where applying a closure by truncation at low P and J can introduce spurious artifacts 13 .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.

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 ∥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. 10 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 29 .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 [46][47][48] .
We present the simulations of a linear device using LAPD physical parameters based on a Helium plasma 29 and a firstof-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 highcollisional 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 ∥ ≃ 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 29,35 and GK 37,38 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 ∼ 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) ∼ (12, 1) in the low collisional regime to (P, J) ∼ (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 simulations [46][47][48] , improvements in our physical model are required to increase the fidelity of our simulations.For example, the inclusion of FLR effects 50,51 is necessary to study small scale fluctuations.A more de- tailed 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 operator 14 .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. 31.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 ∼ 15 eV) are emitted by pulsed plasma discharges in experiments 68 .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 38 .
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.time derivative of T i yields The time derivatives of the GMs contained in Eq. (A1) can be evaluated using the ion GM hierarchy equation given in Eq. ( 12) with (p, j) = (2, 0) and (0, 1).We obtain Using Eq. (A2) with the ion continuity equation given in Eq. (17a) into Eq.(A1), we derive where the collision terms, C 20 i and C 01 i , cancel out due to the energy conservation of the collision operator (see Eq. ( 15)).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. (A3).
We now use the fact that FLR effect are neglected in our model such that N i ≃ n e .Moreover, an expression for the parallel gradient of the parallel ion flow, ∂ z U ∥i , can be obtained from the vorticity equation given in Eq. (22c).At the leadingorder in the long-wavelength limit, one finds that ∂ z J ∥ ≃ 0, such that Finally, introducing the normalized parallel ion heat flux, and using Eq.(A4), the ion temperature equation Eq. (A3) becomes (A6) We remark that Eq. (A6) can also be deduced from other full-F gyrofluid models 45,48 , which consider finite ion temperature effects with a fixed number of GMs.To close Eq.(A6), an expression for the parallel ion heat flux q ∥i must be prescribed.
In contrast to the ion GM hierarchy equation, which treats the ion parallel heat flux q ∥i as a dynamical variable by evolving high-order GMs (see Eq. (A5)), the Braginskii model assumes a high-collisional limit (where λ mp f k ∥ ≪ 1), allowing for the derivation of a collisional closure for q ∥i .Specifically, using the Fokker-Planck collision operator and applying the highcollisional closure via the Chapman-Enskog procedure 52,69,70 , we can obtain the fluid closure q ∥i ≃ −χ ∥ ∂ z T i , where χ ∥ represents the normalized parallel ion thermal conductivity predicted by the Braginskii transport equation 9,70 .Hence, in the high-collisional limit, the ion temperature equation Eq. (A6) reduces to where the temperature anisotropy is neglected such that T ∥i ≃ T i at high-collisionality 12 .Hence, Eq. (A7) 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 16 , 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 71 , resulting in an overestimate of the parallel heat flux, the GM hierarchy equation allows for the evolution of q ∥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.

FIG. 1 .
FIG.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 ≲ 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.

FIG. 2 .
FIG. 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 ∈ [150, 170] (red) and t ∈ [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 ≳ 130.The time t is normalized to R/c s0 .

FIG. 4 .FIG. 5 .
FIG.4.Time-averaged profiles (solid colored lines) of the electrostatic potential ⟨φ ⟩ t,z (left), electron density ⟨n e ⟩ t,z (center), and electron temperature ⟨T e ⟩ 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 ≤ z ≤ 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.

FIG. 9 .
FIG.9.RMS of the electron density n e (normalized to ⟨n e ⟩ t (r)) (top) and electrostatic potential φ (normalized to ⟨φ ⟩ t (r)) (bottom) fluctuations (computed from Fig.7and 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.

FIG. 13 .
FIG. 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 ∥ /v Ti ∈ [1, 2.5] (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).

FIG. 14 .
FIG.14.Absolute value of the GMs N p j (normalized to N 00 ) associated with the distribution functions shown in Fig.13, plotted on a logarithmic scale as a function of p ≥ 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.

FIG. 15 .
FIG. 15.Instantaneous snapshots of the GMs N p j , 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.

2 3 T
∥i U ∥i −U ∥e ∂ z n e n e − ∂ z U ∥e − 2 3 1 n e ∂ z q ∥i + A E n e + (1 − T i ) S N n e .