Gyrokinetic moment-based simulations of the Dimits shift

We present a convergence study of the gyromoment (GM) approach, which is based on projecting the gyrokinetic distribution function onto a Hermite–Laguerre polynomial basis, focused on the cyclone base case (CBC) (Lin et al., Phys. Rev. Lett., vol. 83, no. 18, 1999, pp. 3645–3648) and Dimits shift (Dimits et al., Phys. Plasmas, vol. 7, no. 3, 2000, pp. 969–983) as benchmarks. We report that the GM approach converges more rapidly in capturing the nonlinear dynamics of the CBC than the continuum GENE code (Jenko et al., Phys. Plasmas, vol. 7, no. 5, 2000, pp. 1904–1910) when comparing the number of points representing the velocity space. Increasing the velocity dissipation improves the convergence properties of the GM approach, albeit yielding a slightly larger saturated heat flux. By varying the temperature equilibrium gradient, we show that the GM approach successfully reproduces the Dimits shift (Dimits et al., Phys. Plasmas, vol. 7, no. 3, 2000, pp. 969–983) and effectively captures its width, which is in contrast to the gyrofluid framework. In the collisional regime, the convergence properties of the GM approach improve and a good agreement with previous global particle-in-cell results on transport is obtained (Lin et al., Phys. Rev. Lett., vol. 83, no. 18, 1999, pp. 3645–3648). Finally, we report that the choice of collision model has a minimal impact both on the ion temperature gradient growth rate and on the nonlinear saturated heat flux, at tokamak-relevant collisionality.


Introduction
To reduce the computational cost of gyrokinetic (GK) simulations, gyrofluid (GF) models evolve a limited number of moments of the distribution function, according to conservation laws derived from the GK Boltzmann equation (Grant & Feix 1967;Madsen 2013;Held et al. 2020).However, a comparative study with a set of GK codes presented by Dimits et al. (2000) reveals that the use of GF models can lead to erroneous results.Focusing on the cyclone base case (CBC), a simulation scenario with the parameters of a DIII-D H-mode discharge (Greenfield et al. 1997), and using particle-in-cell (PIC) and continuum GK codes, Dimits et al. (2000) identify a shift between the minimum value of the background temperature gradient that yields a finite saturated nonlinear heat flux, and the linear ion temperature gradient (ITG) stability threshold.In contrast, the GF model of Beer et al. (1995) does not predict the Dimits shift, thus questioning of the validity of GF models for tokamak core simulations.In fact, despite a correct prediction of the ITG linear growth rate at conventional CBC parameters, the GF model does not present a strongly reduced transport in the ITG marginal stability region, overestimating drastically the transport in comparison to GK codes.
In this study, we show that the limitations of the GF models revealed in Dimits et al. (2000) can be overcome with a GK moment-based approach (Frei et al. 2020).This approach extends the drift-kinetic (DK) model presented in Jorge et al. (2017) and is based on the projection of the velocity space dependence of the distribution function on a Hermite-Laguerre polynomial basis at an arbitrary order.This yields an infinite set of fluid-like equations for the basis coefficients, referred to as gyromoment (GMs), that extend the GF models to an arbitrary number of moments.As the number of GK Hermite-Laguerre moments increases, the evolution of the distribution function converges to the one provided by the full-F GK Boltzmann equation.In addition, Jorge et al. (2019) project advanced GK collision operators on the same basis, including the nonlinear GK Fokker-Planck collision operator (Rosenbluth & Longmire 1957), thus providing a GK model that includes an advanced description of collisional effects and combines the accuracy of the GK model with the efficiency of a fluid approach, in particular when a high collisionality regime is considered.Frei et al. (2020) introduce the linear δf flux-tube limit of the GK moment approach, where equilibrium and fluctuating quantities are separated as in the simulation of the CBC and Dimits shift.The δf flux-tube GM model can be considered an extension of δf flux-tube GF models, as the ones introduced by Brizard (1992), Hammett et al. (1992), Beer et al. (1995), Snyder & Hammett (2001) and Scott (2005), thanks to the use of an arbitrary number of GMs to describe the evolution of the perturbations of a Maxwellian equilibrium distribution.Within the δf framework, Frei et al. (2022) also project and use several GK collision operators, in particular the Dougherty model (Dougherty 1964), the Sugama model (Sugama & Watanabe 2006) and the Landau form of the full Coulomb collision model (Rosenbluth et al. 1957;Hazeltine & Meiss 2003).
Using the linear δf GM model, Frei et al. (2023) validate the GM approach with investigations of toroidal ITG modes, trapped electron modes (TEM), microtearing modes (MTM) and kinetic ballooning modes (KBM) in the s − α flux-tube geometry.Hoffmann et al. (2023) present the first nonlinear simulations based on the δf GM model and use them to investigate the evolution of turbulence in a Z-pinch configuration, where the dynamics is dominated by zonal flows (ZF), considering both the collisionless and the collisional cases.The number of GMs needed for convergence in the collisionless case is shown to be less than the number of grid points necessary for convergence in the stateof-the-art continuum GK code GENE (Jenko et al. 2000), depending on the background gradient strength.Moreover, the results reveal that the choice of the collision model significantly impacts the level of saturated transport, even at low collision frequency.
A similar nonlinear δf GM model is presented by Mandell et al. (2018) and Mandell et al. (2022) and implemented in GX, a GPU native code that solves the nonlinear δf GM hierarchy of equations in record times.GX is benchmarked against GS2 (Kotschenreuther et al. 1995), CGYRO (Candy & Waltz 2003), and Stella (Barnes et al. 2019) in linear and nonlinear cases, for tokamak Miller (Miller et al. 1998) and stellarator magnetic equilibria, while considering both adiabatic and kinetic electrons.In contrast to GF models, Mandell et al. (2022) show that the GM approach can retrieve the Dimits shift.
Extending the work in Mandell et al. (2022), the present study aims to investigate the convergence properties of the GM method in the simulation of the CBC and the Dimits shift (Dimits et al. 2000) using the GYACOMO code (Hoffmann et al. 2022).In addition, the impact of the collision model on the CBC is studied by comparing different linear GK collision models, thus extending the work of Lin et al. (1999), which studies heat transport close to the Dimits threshold using a GK particle-in-cell code and a pitchangle collision operator at a physically relevant collision frequency.Finally, the effect of collisions on the Dimits shift is investigated.
Our study confirms that, similarly to GX, the δf GM approach converges to the correct ITG linear stability threshold and nonlinear transport value, which sets it apart from GF models.This is assessed by comparing the results of the GM approach with the Dimits et al. (2000) results as well as results from the GENE code Jenko et al. (2000).We also conduct a detailed convergence study with respect to the velocity space resolution.We analyze the impact of the background temperature gradient levels, the number of GMs, and the level of numerical velocity dissipation on the convergence properties.Finally, we focus on a specific temperature gradient value corresponding to the Dimits regime and investigate the impact of collisions on the level of transport (Lin et al. 1999), as well as the impact of collisions on the Dimits shift.We compare the results obtained using three different GK collision operators: Dougherty (Dougherty 1964), Sugama (Sugama et al. 2009), and Landau (Rosenbluth et al. 1957).
The paper is organized as follows.Section 2 introduces the Hermite-Laguerre GM approach within the context of the CBC.In Sec. 3, we present the benchmarks of the GM results for the ITG threshold and turbulence in the collisionless limit.The convergence study of the Dimits shift is presented in Sec. 4. In Sec. 5, we investigate the effect of collisions.The conclusions follow in Sec. 6.

Gyrokinetic gyromoment approach in a flux-tube configuration
In the present section, we introduce the GK δf model in a field aligned coordinate system.Then, we project the GK Boltzmann equation on a Hermite-Laguerre polynomial basis, thus obtaining a nonlinear gyromoment model, which is an extension of the linear model presented in Frei et al. (2023).Finally, the numerical implementation of the GM model is discussed.

Nonlinear gyrokinetic model
Within the GK approach (Catto 1978;Frieman & Chen 1982;Hazeltine & Meiss 2003), we consider an adiabatic electron model and evolve the ion gyroaveraged distribution function F i (R, v , µ, t), which expresses the probability density of finding an ion with guiding center R, velocity parallel to the magnetic field v , and magnetic moment µ, at a time t.Using the δf approach, we decompose F i as the sum of a time-independent background Maxwellian component and a perturbation g, i.e.
where the Maxwellian distribution is defined as with B(x) = ||B|| the norm of the equilibrium magnetic field, N = N (x) the equilibrium density, T i = T i (x) the equilibrium temperature, m i the ion mass and v 2 thi = 2T i /m i its thermal velocity.The background quantities depend on the particle position x, which is related to the gyrocenter position as x = R + ρ i e ⊥ , where ρ i is the ion Larmor radius and e We assume small fluctuations of the distribution function, g i /F 0 ∼ ∆ ≪ 1, where the scaling parameter ∆ measures the perturbation amplitude relative to the background (Hazeltine & Meiss 2003).We neglect electromagnetic fluctuations and assume that the background electrostatic potential vanishes, therefore denoting with φ(x, t) the perturbed electrostatic potential, such that eφ/T e ∼ ∆ with T e the electron equilibrium temperature.We also neglect MHD equilibrium pressure effects, namely (4π∇P/B 2 )/(b× ∇B/B) ≪ ∆, with P the total pressure.These assumptions yields the nonlinear δf GK Boltzmann equation where we introduce the parallel gradient operator ∇ = b • ∇ and the non-adiabatic part of the ion distribution function perturbation h = g + F 0 q i φ/T i , being q i the ion charge.In Eq. (2.3), φ(R, t) denotes the gyroaveraged electrostatic potential and represents the ion-ion collision term.
2.2.Field aligned coordinate system and magnetic geometry We use the field-aligned coordinates {ψ, α, χ} with ψ the flux surface label, χ the ballooning angle, and α the binormal angle defined as α = q(ψ)χ − ϕ tor , where q(ψ) is the safety factor and ϕ tor the toroidal angle (Hazeltine & Meiss 2003).Within these coordinates, the magnetic field can be expressed in Clebsch form, B = ∇ψ × ∇α.Considering the local limit, we evaluate the equilibrium quantities at the flux tube center position, r 0 , defining B 0 = B(r 0 ) and q 0 = q(ψ(r 0 )).We then normalize the coordinates (Lapillonne et al. 2009;Frei et al. 2023) where R 0 denotes the major radius of the tokamak.The equilibrium magnetic field can thus be written as B = B 0 ∇x × ∇y and the Jacobian of the normalized, field aligned, coordinate system writes where B(z) = B/B 0 is the normalized magnetic field amplitude and b z = ∇z • b.In the field aligned coordinate system, the parallel gradient operator writes (2.6) The perpendicular nonlinear term in Eq. ( 2.3) can be expressed as where f 1 and f 2 are two generic spatially dependent fields and is the Poisson Bracket operator and we introduce the geometric factors being g ij (z) = ∇i • ∇j the metric coefficients for i, j = x, y, z.
In this work we consider the s − α magnetic equilibrium model, which assumes circular concentric flux surfaces using a first order approximation with respect to the inverse aspect ratio ε 0 = r 0 /R 0 (Connor et al. 1978).This includes a constant magnetic shear ŝ, such that the safety factor is expressed as q(x) = q 0 (1 + xŝ).The ballooning angle χ is approximated by the poloidal angle θ, which yields z = θ.Neglecting the MHD equilibrium pressure effects, the metric coefficients are given by   g xx g xy g xz g yx g yy g yz g zx g zy g zz where ε = r 0 (1 + x)/R 0 .The s − α geometry retains finite aspect ratio term in the expression of the magnetic field amplitude B = 1/(1+ε cos z), which yields the derivatives ∂ x ln B = − cos(z) B/R 0 , ∂ y ln B = 0, ∂ z ln B = ǫ sin(z) B, and the Jacobian takes the form J 0 = q 0 B 0 /B.While the s − α model is known to present inconsistencies in the ε ordering (Lapillonne et al. 2009), we consider it here with the purpose of establishing direct comparisons with previous literature results (Lin et al. 1999;Dimits et al. 2000;Frei et al. 2023).

Scale separation
We consider a scale separation between perpendicular and parallel fluctuations, ordering them such that k x ρ s ∼ k y ρ s ∼ k z R 0 , where ρ s = m i c s /(qB 0 ) is the reference sound Larmor radius and c s0 = T e /m i the reference sound speed.Assuming ρ s0 /R 0 ∼ ∆, this implies that k z /k x,y ∼ ∆, which is valid for fluctuations in typical conditions of the core and pedestal of JET and ITER (Giroud et al. 2015).Consequently, we write the magnetic curvature operator as where we neglect terms related to parallel derivatives.The perpendicular nonlinear term present in Eq. ( 2.3) can be simplified by using the same ordering, that is (2.10) Finally, we note that the parallel nonlinear term can be neglected since (2.11) Within these assumptions and considering radially dependent density and temperature equilibrium profiles, the GK Boltzmann equation, Eq. ( 2.3), writes where we introduce the dimensionless velocity coordinates s = v /v thi and w ⊥ = µB/T i , and use the relation In the rest of this work, we express all quantities in dimensionless units according the normalization presented in Table 1, except for the figure labels where we use physical units.In the flux tube limit, we impose constant background gradient profiles, assuming Table 1.Dimensionless variables used in the gyromoment model.For a dimensionless variable A, its equivalent in physical units is explicitly denoted as A ph .
the dimensionless, scale separated, flux tube GK Boltzmann equation writes (2.13) introducing κ N = R 0 /L N the normalized background density gradient, κ T = R 0 /L T the normalized background temperature gradient, τ the ion-electron temperature ratio, q the ion charge number and σ the electron-ion mass ratio.

Nonlinear Hermite-Laguerre-Fourier pseudo spectral formulation
We project the dimensionless GK perturbed distribution function onto a Hermite and Laguerre polynomial basis (Grant & Feix 1967;Madsen 2013;Manas et al. 2017;Adkins & Schekochihin 2018;Staebler et al. 2023) (2.14) where we introduce M pj (x, y, z, t), the ion GM of order (p, j).The Hermite polynomial of order p is defined as The Laguerre polynomial of order j is expressed as and satisfies the orthogonality relation The Laguerre polynomials also satisfy the product identity w ⊥ L j = (2j + 1)L j − (j + 1)L j+1 − jL j−1 (Gradshteyn & Ryzhik 2014).
By using the orthogonality relations, the GMs are obtained from the distribution function as (2.17) The flux tube model considers the ∇x and ∇y directions as periodic, which is valid as long as the domain size is larger than the perpendicular correlation length of the turbulent eddies (Ball et al. 2020).It is thus convenient to express our fields in terms of (k x , k y ) Fourier modes, i.e.
N pj (k x , k y , z, t) = dxdyM pj (x, y, z, t)e −ikxx−ikyy . (2.18) It follows that the Fourier representation of the gyroaveraged electrostatic potential yields with J 0 the Bessel function of the first kind, which can be expressed in terms of Laguerre polynomials as where x +2g xy k x k y +g yy k 2 y (Frei et al. 2020).The projection of Eq. (2.13) onto the Hermite-Laguerre basis yields the gyromoment nonlinear hierarchy in a flux tube configuration, (2.21) In Eq. (2.21), the perpendicular magnetic term, related to the curvature and gradient drifts, writes with C kxky the magnetic curvature operator of Eq. (2.9) expressed in Fourier space, while the parallel magnetic term, related to the Landau damping and the mirror force, is expressed as The background gradient drift terms are for the density and for the temperature.In Eqs.(2.22), (2.23) and (2.25), we also introduce the non-adiabatic ion GM, n pj (k, t) = N pj + qφK j δ p0 /τ .The nonlinear term related to the E × B drift is expressed as where we use the Bessel-Laguerre decomposition, Eq. (2.20), and we express the product of two Laguerre polynomials as a sum of single polynomials using the identity (2.28) Finally, using an adiabatic electron response, we close our system with the dimensionless GK Poisson equation in Fourier space, i.e.
where φ yz is the flux surface average of φ, namely (2.30) While the adiabatic electron model considered here serves as a valuable tool for comparison with previous work, it is important to note its inherent limitations.In fact, evolving the electrons GK equation for the CBC leads to complex small scales phenomena resulting generally, in an increase of linear growth rate and saturated heat flux level due to the electron temperature gradient (ETG) instability (Neiser et al. 2019;Hardman et al. 2022).
To model collisions, we consider the GK Landau operator, which considers the linear limit of exact GK Coulomb collisions, along with the GK Sugama collision model (Sugama et al. 2009) and the GK Dougherty model (Dougherty 1964).The details of these operators and their projection onto the Hermite-Laguerre basis can be found in Frei et al. (2022).For an overview of the qualitative differences between the aforementioned collision operators we also refer to Hoffmann et al. (2023).We set the intensity of the collisions through the nondimensional parameter ν, normalized by the ion-ion collision frequency where ln Λ is the Coulomb Logarithm.

Numerical approach
To solve numerically the GM hierarchy, Eq. (2.21), we evolve a finite set of ion Hermite-Laguerre-Fourier modes, N pj (k x , k y , z, t), with 0 p P and 0 j J.We label a finite set of GMs by the pair (P, J) where P and J represent the maximal polynomial degree of the considered Hermite and Laguerre basis, respectively.For the time integration, we use a standard explicit fourth-order Runge-Kutta time-stepping scheme.

Perpendicular spatial discretization
The Fourier modes are defined for the wavenumbers k x = 2πm/L x , with −N x /2 + 1 m N x /2, and k y = 2πn/L y , with 0 n N y /2 (we exploit the reality condition) being L x and L y the perpendicular dimensions of the flux tube.The nonlinear term, Eq. (2.26), is evaluated using a pseudo-spectral method on each perpendicular plane.More precisely, the spatial derivatives, contained in the Poisson bracket, are obtained in the real space using a backward fast Fourier transform (Frigo & Johnson 2005).The fields are then multiplied in real space and the result transformed back to Fourier space using a forward fast Fourier transform, including a 2/3 anti-aliasing filter (Orszag 1971).To prevent energy pile-up due to the turbulent cascade, we include numerical diffusion of the form µ d (ik/k max ) 4 , where k max is the largest non-aliased allowed wavelength in the simulation and µ d a tunable parameter, set usually to µ d = 0.5.

Parallel spatial discretization
The parallel direction is discretized using a regular grid of N z points with spacing ∆z = 2π/(N z − 1), such that the grid points are z l = l∆z − π for 0 l N z − 1.This constraints our spatial domain to a single poloidal turn around the flux surface (the effect of using a larger number of poloidal turns is described by Ball et al. (2020) and Volčokas et al. (2023)).For evaluating the parallel derivative coming from the Landau damping term in Eq. (2.23), we use a four-point, fourth-order, centered finite differences scheme.
To account for the effect of magnetic shear, we adopt a "twist-and-shift" condition for the parallel boundary conditions, which imposes the coupling between k x and k y modes, i.e.
with A a generic spatially dependent field.For more details on the twist and shift boundary conditions, we refer to the works of Beer et al. (1995) and Ball et al. (2020).
We also add a fourth-order parallel numerical diffusion term of the form µ z /∆z 4 ∂ 4 z , with µ z a tunable parameter, which prevents the decoupling between odd and even points in the z direction (Paruta et al. 2018) and helps to smooth out the oscillations due to the Dirichlet boundary condition applied at the end of the ballooning angle domain.We note that a staggered grid representation can also be implemented to avoid the checkerboard pattern, as it is done for Braginskii's equations in Paruta et al. (2018).

Velocity space discretization and closure
When a finite set of (P, J) GMs is considered, (P + 1) × (J + 1) coupled equations are solved.In order to close the system, the GMs with degree higher than P or J appearing in Eqs.(2.22), (2.23), and (2.26), are assumed to vanish using a truncation closure, i.e.N pj = 0 for p > P or j > J.The truncation closure is the most straightforward to apply and shows good results in Frei et al. (2022), Frei et al. (2023) and Hoffmann et al. (2023).We must note that, in collisionless simulations, the Hermite parallel coupling in Eq. (2.23), combined with the closure by truncation, may lead to recurrence effects, which yields a nonphysical energy transfer from the highest to the lowest Hermite modes.When considering the collisionless limit, we avoid recurrence effects by adding a numerical diffusive term in the velocity space,

Cyclone base case convergence study
We present a benchmark and a convergence analysis of the GM approach in the collisionless CBC (Lin et al. 1999;Dimits et al. 2000).We focus on the ITG linear growth rate and the nonlinear saturated heat flux.We compare our findings with the results obtained from previous work (Dimits et al. 2000), as well as simulations carried out by using the continuum code GENE (Jenko et al. 2000), also considering an adiabatic electron response.We also investigate the convergence of the GM approach by varying the number of evolved moments.
Unless stated otherwise, we employ the conventional CBC parameters: κ T = 6.96, κ N = 2.22, q 0 = 1.4,ǫ = 0.18 and ŝ = 0.8.Regarding the resolution, our linear scan encompasses N x = 8 coupled k x modes and N z = 24 points in the parallel direction.For the nonlinear simulations, we set L x /ρ s = 120, L y /ρ s = 120, and L z /R 0 = 2π, with corresponding grid dimensions N x = 128, N y = 64, and N z = 24 in both the GENE and GM codes.Similarly, a fourth-order parallel dissipation term with µ z = 0.2 is used in both codes (Pueschel et al. 2010).For GENE simulations, we use the standard velocity domain, L v × L µ = 9 × 3, where L v and L µ are the dimensions of the velocity domain along the parallel velocity and magnetic moment directions, respectively.Finally, the velocity fourth-order dissipation parameter of GENE is set to ν v = 0.2 (Pueschel et al. 2010), while the GM simulations use the numerical velocity dissipation frequency η v = 0.001.

Linear convergence and ion temperature gradient stability threshold
To evaluate the ITG growth rate, we numerically solve the moment hierarchy equation, Eq. (2.21), neglecting the nonlinear term in Eq. ( 2.26) and we analyze the time evolution of φ(t) at k x = 0 and z = 0 for various k y values.To capture the parallel dynamics and the coupling due to the magnetic geometry, we evolve 8 k x modes with ∆k x = 2πŝk y .The growth rate γ is determined by fitting the slope of the time evolution of the logarithm of the electrostatic potential amplitude for a given k y mode.
Similarly to our previous work on the entropy mode (Hoffmann et al. 2023), we observe a non-monotonic convergence of the GM growth rates with P and J.The GM approach converges to the correct growth rate with a smaller number of polynomials in the long wavelength limit, k y ρ s ≪ 1, than at short wavelengths.We find that the linear growth rates obtained with the GM approach show good agreement with GENE and Dimits et al. (2000) results when the polynomial basis is sufficiently large, (P, J) (16, 8) (see Fig. 1a).To analyze the convergence properties of the GM approach in more detail, we show the relative error of the peak growth rate (k y ρ s ≃ 0.3), ǫ r = |γ − γ (60,30) |/|γ (60,30) |, in Fig. 1b.The convergence behavior exhibits typical characteristics of spectral methods, that is an exponential and non-monotonic trend with P and J.It is worth noting that Fig. 1b reveals optimal convergence properties when P ∼ 2J, which justifies this choice also in previous works (Frei et al. 2023;Hoffmann et al. 2023;Mandell et al. 2022).This is possibly related to the fact that H n is a Hermite polynomial of degree n in the parallel velocity coordinate, whereas L n is a Laguerre polynomial of degree 2n in the perpendicular velocity coordinate.We note that, while high k y modes converge at a slower rate, they have a reduced impact on the saturated transport level.In fact, according to the mixing length and critical balance estimates (see e.g.Kotschenreuther et al. (1995), Ricci et al. (2006), Schekochihin et al. (2008), Barnes et al. (2011) and Adkins et al. (2023)), the transport of the unstable modes scales as γ/k 2 y .To compare the linear convergence properties between the GM approach and the GENE code, we present a convergence study of the linear growth rate carried out with the GENE code in Fig. 1a.We observe that the GM approach converges faster than the GENE code for small wavenumbers, which can be attributed to the fluid nature of these modes and to a faster convergence of the Bessel-Laguerre decomposition (see Eq. (2.20)).At low velocity resolution, our results show that the GENE code tends to reduce the linear growth rate.This is in contrast to the GM approach, which overestimates the linear growth rates for both the (2, 1) and (4, 2) basis.In addition, Figure 1 presents GENE simulations with a reduced velocity domain, L v × L µ = 4.5 × 1.5 (denoted L v /2).Due to the associated decrease of the velocity grid spacing, we observe an improvement of the results.However, we note that reducing further the velocity domain stabilizes the ITG mode.Similarly, using a N v × N µ = 6 × 4 grid does not yield unstable growth rates when the standard size of the velocity domain is considered.

Nonlinear cyclone base case
We now turn our attention to the nonlinear regime and assess the transport level by examining the normalized ion heat flux Q x .Figure 2a presents the time traces of the heat flux obtained using different GM sets and varying velocity grid resolutions with the GENE code.Both codes converge towards the heat flux value obtained by Dimits et al. (2000).This is also shown in Fig. 2b, where the time-averaged value of the heat flux is shown.We observe that both GENE and the GM codes yield consistent results, exhibiting similar average values and fluctuation amplitudes.With the velocity dissipation frequency η v = 0.001, the GM approach demonstrates faster convergence than the GENE code, providing a good estimate of the heat flux with, approximately, 16 Hermite-Laguerre modes.On the other hand, in agreement with the linear results, we observe that the (4, 2) GM set overestimates the saturated heat flux level.We note that the (2, 1) GM set fails to saturate, leading to a numerical crash.Figure 2b also presents the convergence behavior for two cases with larger numerical velocity dissipation, namely η v = 0.01 and η v = 0.05.In these cases, the moment approach converges more rapidly than with η v = 0.001, but a 30% discrepancy in the converged heat flux value is observed.Since the numerical velocity dissipation term in the GENE code is normalized by the velocity grid size (Pueschel et al. 2010), the incorrect results obtained with the N v × N µ = 8 × 4 and 16 × 8 resolutions is expected to result from the increased strength of the GENE velocity numerical dissipation.However, the simulations where we reduce by a factor two the size of the velocity domain, yield significant changes in the heat flux values (see Fig. 2b) and no clear convergence towards the results obtained with the highest resolutions.Thus, we conclude that the GENE code requires, at least, one hundred grid points in the velocity space to provide a good estimate of the transport level.It is worth noting that, due to recurrence effects, the (P, J) = (2, 1) GM basis presents a pile-up of energy at the low-order Hermite-Laguerre modes, yielding an inaccurate transport level.Alternative closures can potentially prevent this pile up and their use is left for future work.
To further analyze the Hermite-Laguerre representation, we examine the amplitude of the individual GM, defined as E pj = kx,ky,z |N pj (k x , k y , z)|, time-averaged during the quasi-steady saturated state shown in Fig. 3.We observe that the truncation closure used in our simulations affects the Hermite and Laguerre modes differently.In fact, a plateau followed by a sharp decrease can be observed at higher Laguerre degrees j in both the (16, 8) and (30, 15) spectra.On the other hand, a decrease with an approximately constant slope is observed as a function of p.In more details, the (4, 2) GM basis is characterized by a spectrum where the amplitude of the evolved moments is overestimated, although the transport is in good agreement with the converged case.The (8, 4) basis shows a good agreement in the amplitude of the low-degree GMs whereas the (16, 8) simulation shows good agreement with the largest resolution case at all p and j.The reason behing the differences observed between the GM energy spectra in the Hermite and Laguerre directions is twofold.First, the GM hierarchy (2.21) couples the Hermite modes mostly through the magnetic curvature and gradient drifts, connecting the (p, j) GM with the (p ± 2, j) GMs, whereas the Laguerre coupling occurs mostly through neighbouring Laguerre modes, i.e. the (p, j) GM is coupled with the (p, j ± 1) GMs.This explains the oscillations observed in the Hermite modes (p direction in Fig. 3).Second, we recall that the Laguerre polynomials are involved in the Bessel-Laguerre decomposition of Eq.
(2.20), which converges non monotonically with respect to the maximal Laguerre degree, explaining the non monotonic convergence at the highest Laguerre modes in Fig. 3.We note that the convergence of the nonlinear GM simulations is set by the convergence of the peak linear growth rate of the ITG mode, which is in agreement with the findings of our previous work (Hoffmann et al. 2023).In fact, the overestimate of the linear growth rates observed for the (4, 2) basis (see Fig. 1a) corresponds to an increased saturated transport value (see Fig. 2b) and a high amplitude of the GM spectrum on Fig. 3. On the other hand, the (8, 4) basis presents a reduced amplitude for the high degree GMs.

Collisionless Dimits shift
Building upon the study presented in Section 3, we now focus on evaluating the capability of the GM approach to accurately capture the Dimits shift, i.e. the difference between the threshold value of the linear instability and the gradient where a substantial increase of the nonlinear saturated transport level occurs (Dimits et al. 2000).Throughout this section, unless explicitly mentioned, we retain the spatial resolution and numerical dissipation parameters used in Sec. 3.

Collisionless ITG linear threshold
In order to investigate the convergence of the Dimits shift, we define the linear instability threshold as the maximum value of the background temperature gradient for which no growth rates exceeding 0.001 can be observed for 0.05 k y 1.0.This criterion is chosen both due to the inherent difficulty in accurately measuring marginal stability and the acknowledgment that turbulence emerging from such low growth rates, requiring thousands of time units to develop, is expected to have a negligible impact on the plasma dynamics.
Figure 4a presents the results of a parameter scan involving the temperature gradient strength κ T and the number of GMs, having fixed P = 2J, η v = 0.001 and all other parameters as in Sec. 3. Our findings indicate that the GM approach tends to overestimate the ITG threshold when a reduced velocity basis is used, which is in agreement with the observations of Dimits et al. (2000) for GF codes.However, for (P, J) (12, 6), the GM approach yields the threshold value from GK codes reported by Dimits et al. (2000), i.e. κ T ∼ 4. Furthermore, convergence is faster when the gradient level is larger, a trend consistent with Hoffmann et al. (2023).We recall here that this behavior is a consequence of the relative increase of the importance of gradient drift terms (see Eqs. (2.24) and (2.25)), compared to the parallel and perpendicular linear terms (ee Eqs.(2.22) and (2.23)), resulting in reduced coupling of high-degree moments.

Dimits shift and nonlinear convergence study
Turning now to the nonlinear dynamics, we perform a set of simulations to investigate the dependence of the heat transport on the temperature gradient values.We consider the following five sets of GMs: (P, J) = (2, 1), (4, 2), (8, 4), (16, 8), and (30, 15).In addition, we examine the convergence properties of the continuum code GENE by varying the velocity grid resolution, specifically (N v , N µ ) = (8, 4) and (16, 8) with parallel velocity and magnetic moment domain of size L v × L µ = 4.5 × 1.5.We also consider the (N v , N µ ) = (32, 16) and (N v , N µ ) = (42, 24) grids in the velocity domain for L v × L µ = 9 × 3. We measure the radial heat transport by evaluating the ion heat diffusivity, χ = Q x t /(κ T κ N ), with the heat flux averaged over time during the saturated phase of the nonlinear simulations.
Figure 4b displays the heat diffusivity obtained using the GM approach as a function of the background temperature gradient and the number of evolved Hermite-Laguerre moments.By comparing the heat diffusivity with the linear instability threshold determined in the linear case, we also quantify the Dimits shift.In addition, while the critical gradient is influenced by the number of Hermite-Laguerre moments, the width of the shift is not.This provides further support to the observation made in the study of the entropy mode (Hoffmann et al. 2023), which suggests that the constraints on the GM approach convergence are related to the resolution of the primary instability (here the ITG) and not the secondary (the Kelvin-Helmholtz instability) or a possible tertiary one (zonal flows destabilization).
Finally, we investigate the effect of numerical dissipation in velocity space with a convergence study focused on the heat diffusivity, carried out using three different dissipation intensities: η v = 0.05 (Fig. 5a), η v = 0.01 (Fig. 5b) and η v = 0.001 (Fig. 5c).We compare the results with those obtained using GENE (Figure 5d) as well as the results reported in Dimits et al. (2000) for PIC and GK codes.Figure 5 shows that the lowest GM sets, (P, J) = (2, 1) and (4, 2), significantly overestimates the transport values compared Dimits et al. (2000) and GENE, for all dissipation levels and particularly for the lowest background gradient considered.This can be attributed to the overestimate of the linear growth rate, observed in particular at the lowest gradient levels.On the other hand, GENE linear results do not exhibit this overestimate at low resolution, which explains that GENE results do not overestimate the transport at low resolution.
For the (8, 4) and larger GM sets, a good agreement with the results obtained by the various GK codes presented in Dimits et al. (2000) and GENE is observed only for η v = 0.001.In scenarios with low velocity dissipation and close to marginal stability, the GM approach does not outperform GENE in the low resolution limit.This implies the existence of non-Maxwellian velocity space structures in the perturbed distribution function, which are challenging to resolve with only a few Hermite-Laguerre modes.This aligns with our previous study on the Dimits shift in a Z-pinch, where the velocity distributions are compare in more details (Hoffmann et al. 2023).At higher dissipation levels, the GM approach converges faster but yields inaccurate values, in particular when the gradient level becomes smaller than in the CBC.It is worth noting that the nonlinear GM results converge faster as the gradient is increased, similarly to the linear growth rate.

Collisional effects on the Dimits shift
We explore the impact of collisions on the Dimits shift using advanced GK linearized collision operators.We compare the effect of the linear GK Dougherty (DGGK), Sugama (SGGK) and Landau (LDGK) collision operators on the transport level.We first consider the same parameters as the CBC and set the temperature background level to κ T = 5.3.We vary the collision frequency parameter around ν ∼ 0.005, which corresponds to experimental value of the DIII-D discharge considered for the CBC (Lin et al. 1999).We then consider the impact of collisions for different values of κ T , setting η v = 0.
Motivated by the observed correlation between the level of heat transport and the ITG linear growth rate observed in the collisionless case (see Sec. 5), we first study the impact of our different collision models on the linear growth rates.These are presented Fig. 6, the collision ranging from ν = 0.05 to ν = 0.005, and the polynomial basis varying from (P, J) = (4, 2) to (16, 8).In agreement with the collisionless case, we observe that the (4, 2) basis tends to overestimate the growth rate for both collision frequencies across all three collision models.However, faster convergence is observed compared to the collisionless scenario, which is a consequence of the fact that the GK model tends towards a fluid limit at high collisionality.Specifically, the growth rates obtained with the (P, J) = (8, 4) basis closely approach those obtained with the (16, 8) basis.This is not the case in the collisionless regime.Notably, we find that the growth rates do not exhibit significant variations among the different collision operators and are weakly sensitive to the collision frequency for large GM sets.We now perform nonlinear simulations with the GM code at κ T = 5.3 for a set of collision frequencies, namely 0.005 ν 0.05 with the (P, J) = (16, 8) polynomial basis.In Fig. 7a, we present the heat flux time traces obtained by Lin et al. (1999), and the GM approach using the three collision operators considered here.Remarkably, we observe similar results across all collision models, with bursts increasing the transport by a factor five with respect to the baseline level, as also observed by Lin et al. (1999).It is worth mentioning that these bursts are also observed in GK simulations (Kobayashi et al. 2015;Peeters et al. 2016;Hallenbert & Plunk 2022) as well as reduced turbulence modelling (Qi et al. 2020;Ivanov et al. 2020Ivanov et al. , 2022)).Figure 7b shows the average heat diffusivity as a function of collision frequency for the different collision models, compared with the results from Lin et al. (1999).We note a good agreement in terms of trend and transport amplitude,  (1999), reported in gray.The value of the heat diffusivity in the collisionless case is also shown (black dashed line).The error bars represent the standard deviation of the heat flux.
Figure 8.Heat diffusivity as a function of the temperature gradient strength R/LT obtained with the GM approach in the collisionless limit (blue circles) and for ν = 0.005 with the Landau GK collision operator (green diamonds).We report the collisionless results (black dots) of Dimits et al. (2000) and the fit used therein to obtain the threshold value κT = 6 (dashed line).
particularly considering that the results from Lin et al. (1999) are based on global PIC simulations performed with the GTC code (Lin et al. 1998).
Our results show that the choice of collision model does not appear to significantly affect the heat flux at collision frequencies close to the physical values observed in the core of the DIII-D tokamak.The choice of collision models impacts the transport only when the collision frequency is ten times higher than the physically relevant value.This contrasts with Hoffmann et al. (2023), where a noticeable effect of the choice of the collision operator is observed also at low collision frequency.We attribute this difference to two main factors.First, the ITG instability exhibits a stability spectrum with predominantly stable small-scale wavelengths, thus considerably less affected by collisions, with respect to the entropy mode considered by Hoffmann et al. (2023) that develops also on short scales.Second, the collision operators used here account only for ion-ion collisions while electron-ion collisions are considered in the entropy mode investigations.
Finally, we study the effect of the Landau collision operator on the Dimits shift at ν LDGK = 0.005.As shown on Fig. 8, the collisionality tends to smooth out the transition between fully developed turbulence (κ T 6) and zonal flows dominated saturated states (κ T 6), where collisions have a stronger effect on the transport level.This observation is also made in Peeters et al. (2016), where numerical dissipative terms are used to obtain the smoothing of the Dimits shift.This suggests that advanced collision operators are not required to accurately resolve the CBC with an adiabatic electron model, as zonal flow are mostly damped through diffusion in the configuration and velocity spaces.

Conclusions
In the present study, we report on a benchmark and convergence analysis of the GM approach in the local flux-tube δf framework, with a specific focus on the CBC and the Dimits shift.The GM approach converges faster than the GENE code in yielding the correct linear growth rate of the low k y ρ s modes, and therefore a mixing-length estimate of the turbulent transport.In the nonlinear case, our findings highlight that the GM approach accurately captures the nonlinear dynamics of the CBC, while using significantly fewer velocity space points compared to the GENE code.We observe that increasing the intensity of velocity dissipation improves the convergence rate, albeit with a 30% discrepancy in the saturated heat flux value.In addition, the GM approach successfully replicates the GK results reported in Dimits et al. (2000), in contrast to the GF models.The Dimits shift is obtained accurately with a comparable number of moments as the number of GENE velocity grid points, observing slower convergence as the system approaches marginal stability.Nonetheless, it is important to highlight that the GM approach effectively captures the width of the Dimits shift also at low velocity space resolution.This highlights that the primary limitation of the GM approach lies in the convergence of the linear stability threshold rather than in the nonlinear dynamics.Thus, evolving additional GMs in the CBC acts mostly as a fine tuning of the linear ITG instability.It also confirms that a simple model based on the E × B advection, as in Ivanov et al. (2022), is sufficient to predict correctly the dynamics in the Dimits region.Consequently, the balance between Reynolds and diamagnetic stresses, present in the E × B advected temperature assumption (Sarazin et al. 2021), is maintained when increasing the number of evolved GMs in the CBC.However, one must recall that this affirmation is challenged by GK simulations based on the global code GYSELA (Sarazin et al. 2021), suggesting that the disagreement between GYSELA and Ivanov et al. (2022) simulations resides in features of the models such as the local assumption or the boundary conditions.
In the collisional case, our analysis of transport within the Dimits window (κ T = 5.3) reveals a minimal influence of collision on the ITG growth rate around a physical collision frequency of ν ∼ 0.005.This observation holds for all the three GK linear collision models examined, namely Dougherty, Sugama, and Landau.While the GM approach exhibits good agreement with the global PIC results reported in Lin et al. (1999), when varying the collision frequency, the collision model itself does not exhibit a significant impact around the DIII-D core relevant collision frequency.We attribute this result to the longwavelength nature of the ITG instability and to the adiabatic electron model, which prevents the impact of electron-ion collisions on the plasma dynamics.
In a broader context, the present study is a stepping stone towards the efficient modelling of plasma turbulence in the edge region.In fact, as confirmed by the present work, the high collisionality, combined with large background gradients typical of the edge, are ideal conditions to apply the GM approach.

Figure 1 .
Figure 1.(a) CBC linear growth rates, γ obtained by using the GM approach (top) and the GENE code (bottom), compared with the results reported in Dimits et al. (2000) (black crosses).The velocity resolution is scanned by varying the size of the polynomial sets (P, J) for the GM approach and the number of velocity grid points Nv × Nµ for the GENE code.GENE Results with halved velocity domain, L v = 4.5 and Lµ = 1.5, are denoted as Lv/2.(b) Convergence of the relative error ǫr of the CBC linear growth rate obtained with the GM approach at kyρs = 0.3.The error is evaluated with respect to the growth rate evaluated with (P, J) = (60, 30) GM, γ (60,30) , namely ǫr = |γ − γ (60,30) |/|γ (60,30) | .

Figure 2 .
Figure2.(a) Time traces of the heat flux for the CBC with comparison between the gyromoment approach, GENE and the result ofDimits et al. (2000), Qx ≈ 35.Various velocity resolutions are used, in particular different GM sets (P, J) and numbers of points in the GENE velocity grid Nv × Nµ.(b) Convergence in the saturated heat flux value with respect to the number of points representing the velocity space N v p = (P +1)×(J +1) for the GM approach and N v p = Nv ×Nµ for GENE.The configuration space resolution is Nx = 128, Ny = 64 and Nz = 24 for both the GENE code and the GM approach.The error bars reflect the average fluctuation amplitude around the time-averaged transport value.The GENE simulations denoted by Lv/2 are obtained with a halved velocity domain size.

Figure 4 .
Figure 4. (a) ITG growth rate at kyρs = 0.3 for different background temperature value R0/LT and (P, J), with J = P/2, R0/LN = 2.22 and ηv = 0.001.(b) Heat diffusivity obtained with the GM approach and the threshold for the ITG linear growth rate (black dashed line).The gap in the R0/LT values between the ITG linear threshold and the non-zero χ values represents the Dimits shift.

Figure 7 .
Figure7.Study of collision effects on the heat flux with the GM approach and the Dougherty (orange), Sugama (blue) and Landau (green) operators for κT = 5.3 for ν = 0.005.The heat flux time traces (a) and the average heat diffusivity (b) are compared with the results fromLin et al. (1999), reported in gray.The value of the heat diffusivity in the collisionless case is also shown (black dashed line).The error bars represent the standard deviation of the heat flux.
Shukla et al. (2022)6)2011)nts is still an open issue, more sophisticated closure models such as the semi-collisional closure proposed byZocco & Schekochihin (2011),Loureiro et al. (2016)exist and might reduce the appearance of recurrence effects without introducing an artificial dissipation term.One can also mention the work ofShukla et al. (2022)which presents a method for formulating closures that learn from kinetic simulation data.
2.33)with η v a tunable parameter.We note that this term, similar to a Dougherty collision operator, has the advantage of conserving mass, momentum, and energy.In this work, we ensure that this artificial term does not impact our results by comparing them at different values of η v (see Sec. 4.2).While a generalization of the Braginskii's equation closure to