Gyrokinetic moment-based simulations of the Dimits shift

Abstract

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.

1. 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, Wiesenberger & Kendl 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, Cowley & Hammett (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 with 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, Jorge & Ricci 2020). This approach extends the drift-kinetic model presented in Jorge, Ricci & Loureiro (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 gyromoments (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, Frei & Ricci (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 $\delta 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 $\delta f$ flux-tube GM model can be considered an extension of $\delta f$ flux-tube GF models, such as the ones introduced by Brizard (1992), Hammett, Dorland & Perkins (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 $\delta f$ framework, Frei, Hoffmann & Ricci (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, MacDonald & Judd 1957; Hazeltine & Meiss 2003).

Using the linear $\delta f$ GM model, Frei et al. (2023) validate the GM approach with investigations of toroidal ITG modes, trapped electron modes, microtearing modes and kinetic ballooning modes in the $s-\alpha$ flux-tube geometry. Hoffmann, Frei & Ricci (2023) present the first nonlinear simulations based on the $\delta 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, 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 state-of-the-art continuum GK code GENE (Jenko, Dorland & Kotschenreuther 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 $\delta f$ GM model is presented by Mandell, Dorland & Landreman (2018) and Mandell et al. (2022) and implemented in GX, a GPU native code that solves the nonlinear $\delta f$ GM hierarchy of equations in record times. The GX code is benchmarked against GS2 (Kotschenreuther, Rewoldt & Tang 1995), CGYRO (Candy & Waltz 2003) and Stella (Barnes, Parra & Landreman 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 open source MPI-based GYACOMO code (Hoffmann, Frei & Ricci 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 PIC code and a pitch-angle 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 $\delta 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 analyse 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 (1964), Sugama, Watanabe & Nunami (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 § 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 § 4. In § 5, we investigate the effect of collisions. The conclusions follow in § 6.

2. Gyrokinetic gyromoment approach in a flux-tube configuration

In the present section, we introduce the GK $\delta 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 GM 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.

2.1. 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 $\langle F_i\rangle (\boldsymbol R, v_\parallel, \mu, t)$, which expresses the probability density of finding an ion with guiding centre $\boldsymbol R$, velocity parallel to the magnetic field $v_\parallel$ and magnetic moment $\mu$, at a time $t$. Using the $\delta f$ approach, we decompose $\langle F_i\rangle$ as the sum of a time-independent background Maxwellian component and a perturbation $g$, i.e.

(2.1)\begin{equation} \langle F_i\rangle(\boldsymbol R, v_\parallel, \mu, t) = F_0(\boldsymbol x, v_\parallel, \mu) + g(\boldsymbol R, v_\parallel, \mu, t), \end{equation}

where the Maxwellian distribution is defined as

(2.2)\begin{equation} F_0(\boldsymbol x, v_\parallel, \mu) =\frac{N}{({\rm \pi}^{1/2}v_{\rm thi})^{3}} \exp \left(-\frac{m_i v_\parallel^2}{2T_i}-\frac{\mu B}{T_i}\right), \end{equation}

with $B(\boldsymbol x)=\|\boldsymbol B\|$ the norm of the equilibrium magnetic field, $N=N(\boldsymbol x)$ the equilibrium density, $T_i=T_i(\boldsymbol x)$ the equilibrium temperature, $m_i$ the ion mass and $v_{\rm thi}^2 = 2T_i/m_i$ its thermal velocity. The background quantities depend on the particle position $\boldsymbol x$, which is related to the gyrocentre position as $\boldsymbol x = \boldsymbol R + \rho _i \boldsymbol e_\perp$, where $\rho _i$ is the ion Larmor radius and $\boldsymbol e_\perp =[\boldsymbol R - (\boldsymbol R \boldsymbol {\cdot } \boldsymbol b)\boldsymbol {b}] /\|\boldsymbol R - (\boldsymbol R \boldsymbol {\cdot } \boldsymbol b)\boldsymbol {b}\|$, with $\boldsymbol {b}=\boldsymbol B/B$.

We assume small fluctuations of the distribution function, $g_i/F_0\sim \varDelta \ll 1$, where the scaling parameter $\varDelta$ 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 $\phi (\boldsymbol x,t)$ the perturbed electrostatic potential, such that $e\phi /T_e\sim \varDelta$, with $T_e$ the electron equilibrium temperature. We also neglect magnetohydrodynamic (MHD) equilibrium pressure effects, namely $(4{\rm \pi} \boldsymbol {\nabla } P/B^2)/(\boldsymbol {b}\times \boldsymbol {\nabla } B/B)\ll \varDelta$, with $P$ the total pressure. These assumptions yield the nonlinear $\delta f$ GK Boltzmann equation

(2.3)\begin{align} & \partial_t g - \frac{q_i}{m_i}\nabla_\parallel\bar \phi\partial_{v\parallel} g + \frac{1}{B}(\boldsymbol{b}\times\boldsymbol{\nabla}\bar \phi)\boldsymbol{\cdot} \boldsymbol{\nabla} g +\frac{m_i}{q_i B} \left(v_\parallel^2 + \frac{\mu B}{m_i}\right)(\boldsymbol{b}\times\boldsymbol{\nabla} \ln B)\boldsymbol{\cdot} \boldsymbol{\nabla} h +v_\parallel\nabla_\parallel h\nonumber\\ & \quad - \frac{\mu B}{m_i}\nabla_\parallel \ln B\partial_{v\parallel} h +\frac{1}{B}\boldsymbol{b}\times\left(\frac{\boldsymbol{\nabla} N}{N} +\left[ \frac{m_i v_\parallel^2}{2T_i}+\frac{\mu B}{T_i}-\frac{3}{2} \right]\frac{\boldsymbol{\nabla} T_i}{T_i}\right)\boldsymbol{\cdot} \boldsymbol{\nabla}\bar \phi = C_{\rm ii}, \end{align}

where we introduce the parallel gradient operator $\nabla _\parallel = \boldsymbol {b}\boldsymbol {\cdot } \boldsymbol {\nabla }$ and the non-adiabatic part of the ion distribution function perturbation $h=g+ F_0 q_i\phi /T_i$, with $q_i$ the ion charge. In (2.3), $\bar \phi (\boldsymbol R,t)$ denotes the gyroaveraged electrostatic potential and $C_{\rm ii} = C_{\rm ii}(\boldsymbol R, v_\parallel, \mu )$ represents the ion–ion collision term.

2.2. Field-aligned coordinate system and magnetic geometry

We use the field-aligned coordinates $\{\psi,\alpha,\chi \}$, with $\psi$ the flux surface label, $\chi$ the ballooning angle and $\alpha$ the binormal angle defined as $\alpha = q(\psi )\chi - \varphi _{\rm tor}$, where $q(\psi )$ is the safety factor and $\varphi _{\rm tor}$ the toroidal angle (Hazeltine & Meiss 2003). Within these coordinates, the magnetic field can be expressed in Clebsch form, $\boldsymbol B = \nabla \psi \times \nabla \alpha$. Considering the local limit, we evaluate the equilibrium quantities at the flux-tube centre position, $r_0$, defining $B_0=B(r_0)$ and $q_0 = q(\psi (r_0))$. We then normalize the coordinates (Lapillonne et al. 2009; Frei et al. 2023)

(2.4a–c)\begin{equation} x=\frac{q_0}{r_0 B_0}[\psi-\psi(r_0)],\quad y=\frac{R_0}{q_0}\alpha,\quad z=\chi, \end{equation}

where $R_0$ denotes the major radius of the tokamak. The equilibrium magnetic field can thus be written as $\boldsymbol B = B_0\boldsymbol {\nabla } x\times \boldsymbol {\nabla } y$ and the Jacobian of the normalized field-aligned coordinate system writes

(2.5)\begin{equation} J_{xyz} = (\boldsymbol{\nabla} z \boldsymbol{\cdot} \boldsymbol{\nabla} x \times \boldsymbol{\nabla} y)^{{-}1} = B_0/(\boldsymbol{\nabla} z \boldsymbol{\cdot} \boldsymbol B) = (\hat B b_z)^{{-}1}, \end{equation}

where $\hat B(z) = B/B_0$ is the normalized magnetic field amplitude and $b_z= \boldsymbol {\nabla } z \boldsymbol {\cdot } \boldsymbol {b}$. In the field-aligned coordinate system, the parallel gradient operator writes

(2.6)\begin{equation} \nabla_\parallel{=} \frac{1}{J_{xyz}\hat B}\partial_z. \end{equation}

The perpendicular nonlinear term in (2.3) can be expressed as

(2.7)\begin{equation} (\boldsymbol{b}\times \boldsymbol{\nabla} f_1)\boldsymbol{\cdot} \boldsymbol{\nabla} f_2 = \hat B^{{-}1}[\varGamma_1 \{f_1,f_2\}_{xy} + \varGamma_2 \{f_1,f_2\}_{xz} + \varGamma_3 \{f_1,f_2\}_{yz}], \end{equation}

where $f_1$ and $f_2$ are two generic spatially dependent fields and $\{f_1,f_2\}_{ij}=\partial _i f_1 \partial _j f_2 - \partial _j f_1 \partial _i f_2$ is the Poisson bracket operator and we introduce the geometric factors

(2.8)\begin{equation} \left.\begin{array}{c} \varGamma_1=g^{xx}g^{yy} - g^{xy}g^{xy},\\ \varGamma_2=g^{xx}g^{yz} - g^{xy}g^{xz},\\ \varGamma_3=g^{xy}g^{yz} - g^{yy}g^{xz}, \end{array}\right\} \end{equation}

with $g^{ij}(z) = \boldsymbol {\nabla } i \boldsymbol {\cdot } \boldsymbol {\nabla } j$ the metric coefficients for $i,j=x,y,z$.

In this work, we consider the $s-\alpha$ magnetic equilibrium model, which assumes circular concentric flux surfaces using a first-order approximation with respect to the inverse aspect ratio $\varepsilon _0=r_0/R_0$ (Connor, Hastie & Taylor 1978). This includes a constant magnetic shear $\hat s$, such that the safety factor is expressed as $q(x)=q_0(1+x\hat s)$. The ballooning angle $\chi$ is approximated by the poloidal angle $\theta$, which yields $z=\theta$. Neglecting the MHD equilibrium pressure effects, the metric coefficients are given by

(2.9)\begin{equation} \begin{pmatrix} g_{xx} & g_{xy} & g_{xz}\\ g_{yx} & g_{yy} & g_{yz}\\ g_{zx} & g_{zy} & g_{zz} \end{pmatrix} = \begin{pmatrix} 1 & \hat sz & 0\\ \hat sz & 1+(\hat sz)^2 & \varepsilon^{{-}1}\\ 0 & \varepsilon^{{-}1} & \varepsilon^{{-}2} \end{pmatrix}, \end{equation}

where $\varepsilon =r_0(1+x)/R_0$. The $s-\alpha$ geometry retains the finite aspect ratio term in the expression of the magnetic field amplitude $\hat B = 1/(1+\varepsilon \cos z)$, which yields the derivatives $\partial _x \ln B = -\cos (z)\hat B/R_0$, $\partial _y \ln B = 0$, $\partial _z \ln B = \epsilon \sin (z)\hat B$ and the Jacobian takes the form $J_0=q_0 B_0/B$. While the $s-\alpha$ model is known to present inconsistencies in the $\varepsilon$ 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).

2.3. Scale separation

We consider a scale separation between perpendicular and parallel fluctuations, ordering them such that $k_{x}\rho _{s} \sim k_{y}\rho _{s} \sim k_z R_0$, where $\rho _{s} = m_i c_{s}/(q B_0)$ is the reference sound Larmor radius and $c_{s0}=\sqrt {T_e/m_i}$ the reference sound speed. Assuming $\rho _{s0}/R_0 \sim \varDelta$, this implies that $k_z/k_{x,y}\sim \varDelta$, 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 $(\boldsymbol {b}\times \boldsymbol {\nabla } \ln B)\boldsymbol {\cdot } \boldsymbol {\nabla } = \varGamma _1\hat B^{-1} \mathcal {C}_{xy}$, with

(2.10)\begin{equation} \mathcal{C}_{xy} ={-}\left[ \partial_y\ln B + \frac{\varGamma_2}{\varGamma_1}\partial_z\ln B \right] \frac{\partial}{\partial x} + \left[ \partial_x\ln B - \frac{\varGamma_3}{\varGamma_1}\partial_z\ln B \right] \frac{\partial}{\partial y}, \end{equation}

where we neglect terms related to parallel derivatives. The perpendicular nonlinear term present in (2.3) can be simplified by using the same ordering, that is

(2.11)\begin{equation} (\boldsymbol{b}\times\boldsymbol{\nabla}\bar \phi)\boldsymbol{\cdot} \boldsymbol{\nabla} g = \varGamma_1\hat B^{{-}1}\{\bar \phi,g\}_{xy}. \end{equation}

Finally, we note that the parallel nonlinear term can be neglected since

(2.12)\begin{equation} (q/m) b_z \partial_z \bar \phi \partial_{v\parallel} g \sim \bar \phi g/(B_0R_0\rho_{s0})\ll \varDelta. \end{equation}

Within these assumptions, and considering radially dependent density and temperature equilibrium profiles, the GK Boltzmann equation, (2.3), writes

(2.13)\begin{align} & \partial_t g + \frac{\hat B}{B} \{\bar \phi,g\}_{xy}+ \frac{\hat B}{B}\frac{m_i v_{thi}^2}{2q_i}\ (2s_\parallel^2+w_\perp) \mathcal{C}_{xy}h \nonumber\\ & \quad + \frac{\hat B}{B} \left[ \partial_x\ln N + \left(s_\parallel^2+w_\perp{-}\frac{3}{2}\right)\partial_x \ln T_{i} \right]F_0\partial_y\bar \phi\nonumber\\ & \quad + \frac{v_{\rm thi}}{2J_{xyz}\hat B}[2s_\parallel \partial_z h - w_\perp \partial_z\ln B\partial_{s\parallel} h] = C_{\rm ii}, \end{align}

where we introduce the dimensionless velocity coordinates $s_\parallel = v_\parallel /v_\textrm {thi}$ and $w_\perp =\mu B/T_i$, and use the relation $\varGamma _1=\hat B^2$.

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 $A(\boldsymbol x)\sim A(0)=A_0$ and $|\boldsymbol {\nabla } A| \sim |A_0/L_A|$ for a generic background quantity $A$. Finally, the dimensionless, scale separated, flux-tube GK Boltzmann equation writes

(2.14)\begin{align} & \partial_t g+ \{\bar \phi,g\}_{xy}+ \frac{\tau }{q}(2 s_\parallel^2+w_\perp ) \mathcal{C}_{xy} h +\left[ \kappa_N + \left(s_\parallel^2+w_\perp{-}\frac{3}{2}\right) \kappa_T \right]\partial_y \bar \phi\nonumber\\ & \quad + \frac{\sqrt{2}}{2}\sqrt{\tau}\frac{1}{J_{xyz}\hat B} [2s_\parallel \partial_z h - w_\perp\partial_z\ln B\partial_{s\parallel} h] = C_{\rm ii}, \end{align}

introducing $\kappa _N={R_0}/{L_N}$ the normalized background density gradient, $\kappa _T={R_0}/{L_T}$ the normalized background temperature gradient, $\tau$ the ion–electron temperature ratio, $q$ the ion charge number and $\sigma$ the electron–ion mass ratio.

Table 1. Dimensionless variables used in the GM model. For a dimensionless variable $A$, its equivalent in physical units is explicitly denoted as $A^\textrm {ph}$.

2.4. Nonlinear Hermite–Laguerre–Fourier pseudospectral 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, Belli & Candy 2023)

(2.15)\begin{equation} g(x,y,z,s_\parallel,w_\perp,t) = \sum_{p=0}^\infty\sum_{j=0}^\infty M^{pj}(x,y,z,t) H_p(s_\parallel)L_j(w_\perp), \end{equation}

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

(2.16)\begin{equation} H_p(s_\parallel)=\frac{({-}1)^p }{\sqrt{2^p p!}}\, {\rm e}^{s_\parallel^2} \frac{{\rm d}^p}{ {\rm d}s_\parallel^p} \,{\rm e}^{{-}s_\parallel^2}, \end{equation}

which is a normalized version of the physicist's Hermite polynomials such that $\int _{-\infty }^\infty \mathrm d s_\parallel H_p H_{p'} \,\textrm {e}^{-s_\parallel ^2}=\sqrt {{\rm \pi} }\delta _{pp'}$, where $\delta _{pp'}$ denotes the Kronecker delta. With this definition, useful Hermite product and derivation identities follow, such as $s_\parallel H_p = \sqrt {(p+1)/2}H_{p+1} + \sqrt {p/2}H_{p-1}$ and $\partial _{s\parallel } H_p =\sqrt {2p}H_{p-1}$. The Laguerre polynomial of order $j$ is expressed as

(2.17)\begin{equation} L_j(w_\perp)=\frac{{\rm e}^{w_\perp}}{j!}\frac{{\rm d}^j}{{\rm d}w_\perp^j}w_\perp^j \,{\rm e}^{{-}w_\perp} ,\end{equation}

and satisfies the orthogonality relation $\int _{0}^\infty \mathrm d w_\perp L_j L_{j'} \,\textrm {e}^{-w_\perp }=\delta _{jj'}$. The Laguerre polynomials also satisfy the product identity $w_\perp 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.18)\begin{equation} M^{pj}(x,y,z,t) = \iint \mathrm d w_\perp \,\mathrm d s_\parallel g(x,y,z,s_\parallel,w_\perp,t) H_p(s_\parallel) L_j(w_\perp). \end{equation}

The flux-tube model considers the $\boldsymbol {\nabla } x$ and $\boldsymbol {\nabla } 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, Brunner & Ajay 2020). It is thus convenient to express our fields in terms of $(k_x,k_y)$ Fourier modes, i.e.

(2.19)\begin{equation} N^{pj}(k_x,k_y,z,t) = \iint \mathrm{d} x \,\mathrm{d} y M^{pj}(x,y,z,t) \,{\rm e}^{-{\rm i} k_x x-{\rm i} k_y y }. \end{equation}

It follows that the Fourier representation of the gyroaveraged electrostatic potential yields

(2.20)\begin{equation} J_0(\sqrt{l w_\perp})\phi(k_x,k_y,z,t) = \iint \mathrm{d} x \,\mathrm{d} y \bar \phi(\boldsymbol x,t)\, {\rm e}^{-{\rm i} k_x x-{\rm i} k_y y }, \end{equation}

with $J_0$ the Bessel function of the first kind, which can be expressed in terms of Laguerre polynomials as

(2.21)\begin{equation} J_0(\sqrt{l w_\perp}) = \sum_{n=0}^\infty \mathcal{K}_n(l)L_n(w_\perp),\end{equation}

where $\mathcal {K}_n(l)=l^n \,\textrm {e}^{-l}/n!$, $l(k_x,k_y,z)=\tau k_\perp ^2/2$ and the perpendicular wavenumber $k_\perp ^2 = g^{xx}k_x^2 + 2 g^{xy}k_x k_y + g^{yy} k_y^2$ (Frei et al. 2020). The projection of (2.14) onto the Hermite–Laguerre basis yields the GM nonlinear hierarchy in a flux-tube configuration

(2.22)\begin{equation} \partial_t N^{pj} + \mathcal{S}^{pj} + \mathcal{M}_{{\perp}}^{pj} + \mathcal{M}_{{\parallel}}^{pj} + \mathcal{D}_{N}^{pj} + \mathcal{D}_{T}^{pj} = \mathcal{C}_{\rm ii}^{pj}. \end{equation}

In (2.22), the perpendicular magnetic term, related to the curvature and gradient drifts, writes

(2.23)\begin{align} \mathcal{M}_{{\perp}}^{pj} & = \frac{\tau}{q} \mathcal{C}_{k_x k_y} [\sqrt{(p+1)(p+2)} n^{p+2,j} + (2p+1)n^{pj} + \sqrt{p(p-1)}n^{p-2,j}] \nonumber\\ & \quad + \frac{\tau}{q} \mathcal{C}_{k_x k_y} [(2j+1)n^{pj} - (j+1)n^{p,j+1}-jn^{p,j-1}], \end{align}

with $C_{k_x k_y}$ the magnetic curvature operator of (2.10) expressed in Fourier space, while the parallel magnetic term, related to the Landau damping and the mirror force, is expressed as

(2.24)\begin{align} \mathcal{M}_{{\parallel}}^{pj} & = \frac{\hat B^{{-}1}}{J_{xyz}} \sqrt{\tau}\{\partial_z \aleph^{p\pm1,j} - \partial_z\ln B [(j+1)\aleph^{p\pm1,j}-j\aleph^{p\pm1,j-1}] \nonumber\\ & \quad+\partial_z\ln B\sqrt{p}[(2j+1)n^{p-1,j} -(j+1)n^{p-1,j+1} - jn^{p-1,j-1}]\}, \end{align}

with $\aleph ^{p\pm 1,j}=\sqrt {p+1} n^{p+1,j} + \sqrt {p} n^{p-1,j}$. The background gradient drift terms are

(2.25)\begin{equation} \mathcal{D}_{N}^{pj} = ik_y\kappa_N\mathcal{K}_j\phi\delta_{p0},\end{equation}

for the density and

(2.26)\begin{equation} \mathcal{D}_{T}^{pj} = ik_y\phi\kappa_T \left\{\mathcal{K}_j\left[ \frac{1}{\sqrt{2}}\delta_{p2} -\delta_{p0} \right]+ [(2j+1)\mathcal{K}_j-(j+1)\mathcal{K}_{j+1}-j\mathcal{K}_{j-1}]\delta_{p0}\right\},\end{equation}

for the temperature. In (2.23), (2.24) and (2.26), we also introduce the non-adiabatic ion GM, $n^{pj}(\boldsymbol k,t)=N^{pj}+q \phi \mathcal {K}_j\delta _{p0}/\tau$.

The nonlinear term related to the $\boldsymbol E\times \boldsymbol B$ drift is expressed as

(2.27)\begin{equation} \mathcal{S}^{pj} = \sum_{n=0}^{\infty}\left\{\mathcal{K}_n\phi,\sum_{s=0}^{n+j}d_{njs} N^{ps} \right\}_{xy}, \end{equation}

where we use the Bessel–Laguerre decomposition, (2.21), and we express the product of two Laguerre polynomials as a sum of single polynomials using the identity

(2.28)\begin{equation} L_jL_n=\sum_{s=0}^{n+j}d_{njs}L_s,\end{equation}

with

(2.29)\begin{equation} d_{njs} = \sum_{n_1=0}^n\sum_{j_1=0}^j\sum_{s_1=0}^s \frac{({-}1)^{n_1+j_1+s_1}}{n_1!j_1!s_1!}\binom{n}{n_1}\binom{j}{j_1}\binom{s}{s_1}. \end{equation}

Finally, using an adiabatic electron response, we close our system with the dimensionless GK Poisson equation in Fourier space, i.e.

(2.30)\begin{equation} \left[ 1 + \frac{q^2}{\tau}\left(1-\sum_{n=0}^{\infty}\mathcal{K}^2_n\right)\right]\phi - \langle \phi \rangle_{yz} = q\sum_{n=0}^{\infty}\mathcal{K}_n N^{0n}, \end{equation}

where $\langle \phi \rangle _{yz}$ is the flux surface average of $\phi$, namely

(2.31)\begin{equation} \langle \phi \rangle_{yz} = \frac{1}{\displaystyle \int\mathrm dz J_{xyz}}\int\mathrm dz J_{xyz}\phi(k_x,k_y,z,t)\delta_{k_y0}. \end{equation}

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-scale phenomena resulting, generally, in an increase of linear growth rate and saturated heat flux level due to the electron temperature gradient 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. (2021). 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 non-dimensional parameter $\nu$, normalized by the ion–ion collision frequency

(2.32)\begin{equation} \nu_{\rm ii} = \frac{4\sqrt{\rm \pi}}{3}\frac{R_0 N e^4 \ln \varLambda}{c_s m_i^{1/2}T_i^{3/2}}, \end{equation}

where $\ln \varLambda$ is the Coulomb Logarithm.

2.5. Numerical approach

To solve numerically the GM hierarchy, (2.22), we evolve a finite set of ion Hermite–Laguerre–Fourier modes, $N^{pj}(k_x,k_y,z,t)$, with $0\leq p \leq P$ and $0\leq j \leq 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.

2.5.1. Perpendicular spatial discretization

The Fourier modes are defined for the wavenumbers $k_x = 2{\rm \pi} m/L_x$, with $-N_{x}/2 +1\leq m \leq N_{x}/2$, and $k_y=2{\rm \pi} n/L_y$, with $0\leq n \leq N_{y}/2$ (we exploit the reality condition), where $L_x$ and $L_y$ are the perpendicular dimensions of the flux tube. The nonlinear term, (2.27), 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 $\mu _d(ik/k_{\max })^4$, where $k_{\max }$ is the largest non-aliased allowed wavelength in the simulation and $\mu _d$ a tuneable parameter, set usually to $\mu _d =0.5$.

2.5.2. Parallel spatial discretization

The parallel direction is discretized using a regular grid of $N_z$ points with spacing ${\rm \Delta} z=2{\rm \pi} /(N_z-1)$, such that the grid points are $z_l = l {\rm \Delta} z - {\rm \pi}$ for $0\leq l \leq N_z-1$. This constrains 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; Volčokas, Ball & Brunner 2023). For evaluating the parallel derivative coming from the Landau damping term in (2.24), we use a four-point fourth-order centred 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.

(2.33)\begin{equation} A(k_x,k_y,z) = A(k_x+2{\rm \pi} \hat s k_y,k_y,z+2{\rm \pi}), \end{equation}

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 $\mu _z/{\rm \Delta} z^4 \partial _z^4$, with $\mu _z$ a tuneable 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 is done for Braginskii's equations in Paruta et al. (2018).

2.5.3. Velocity space discretization and closure

When a finite set of $(P,J)$ GMs is considered, $(P+1)\times (J+1)$ coupled equations are solved. In order to close the system, the GMs with degree higher than $P$ or $J$ appearing in (2.23), (2.24) and (2.27), 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 (2.24), combined with the closure by truncation, may lead to recurrence effects, which yields a non-physical 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

(2.34)\begin{align} D_v^{pj} & ={-}\eta_{v} (p + 2j) N^{pj} +\eta_{v} [ N^{10}\delta_{p1}\delta_{j0} - \tfrac{2}{3} (\sqrt{2}N^{01}-N^{20})\delta_{p2}\delta_{j0} \nonumber\\ & \quad - \tfrac{2}{3} (\sqrt{2}N^{20}-N^{01})\delta_{p0}\delta_{j1}], \end{align}

with $\eta _v$ a tuneable 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 $\eta _{v}$ (see § 4.2). While a generalization of the Braginskii equation closure to an arbitrary number of moments is still an open issue, more sophisticated closure models such as the semi-collisional closure proposed by Zocco & Schekochihin (2011) and 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 of Shukla et al. (2022), which presents a method for formulating closures that learn from kinetic simulation data.

3. 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: $\kappa _T=6.96$, $\kappa _N=2.22$, $q_0=1.4$, $\epsilon =0.18$ and $\hat s = 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/\rho _s=120$, $L_y/\rho _s=120$ and $L_z/R_0=2{\rm \pi}$, 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 $\mu _z=0.2$ is used in both codes (Pueschel, Dannert & Jenko 2010). For GENE simulations, we use the standard velocity domain, $L_{v_\parallel }\times L_\mu = 9\times 3$, where $L_{v_\parallel }$ and $L_{\mu }$ 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 $\nu _v=0.2$ (Pueschel et al. 2010), while the GM simulations use the numerical velocity dissipation frequency $\eta _{v}=0.001$.

3.1. Linear convergence and ion temperature gradient stability threshold

To evaluate the ITG growth rate, we numerically solve the moment hierarchy equation, (2.22), neglecting the nonlinear term in (2.27), and we analyse the time evolution of $\phi (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 ${\rm \Delta} k_x = 2{\rm \pi} \hat s k_y$. The growth rate $\gamma$ 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\rho _s\ll 1$, than at short wavelengths. We find that the linear growth rates obtained with the GM approach show good agreement with the GENE and Dimits et al. (2000) results when the polynomial basis is sufficiently large, $(P,J)\gtrsim (16,8)$ (see figure 1a).

Figure 1. (a) CBC linear growth rates, $\gamma$ 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 $N_{v_\parallel }\times N_\mu$ for the GENE code. GENE results with halved velocity domain, $L_{v\parallel }=4.5$ and $L_\mu =1.5$, are denoted as $L_v/2$. (b) Convergence of the relative error $\epsilon _r$ of the CBC linear growth rate obtained with the GM approach at $k_y\rho _s=0.3$. The error is evaluated with respect to the growth rate evaluated with $(P,J)=(60,30)$ GM, $\gamma _{(60,30)}$, namely $\epsilon _r=|\gamma -\gamma _{(60,30)}|/|\gamma _{(60,30)}|$ .

To analyse the convergence properties of the GM approach in more detail, we show the relative error of the peak growth rate ($k_y\rho _s\simeq 0.3$), $\epsilon _r=|\gamma -\gamma _{(60,30)}|/|\gamma _{(60,30)}|$, in figure 1(b). The convergence behaviour exhibits typical characteristics of spectral methods, that is an exponential and non-monotonic trend with $P$ and $J$. It is worth noting that figure 1(b) reveals optimal convergence properties when $P\sim 2J$, which justifies this choice also in previous works (Mandell et al. 2022; Frei et al. 2023; Hoffmann et al. 2023). 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, Parra & Schekochihin 2011; Adkins, Ivanov & Schekochihin 2023), the transport of the unstable modes scales as $\gamma /k_y^2$.

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 figure 1(a). 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 (2.21)). At low velocity resolution, our results show that the GENE code tends to reduce the linear growth rates. This is in contrast to the GM approach, which overestimates the linear growth rates for both the $(2,1)$ and $(4,2)$ bases. In addition, figure 1 presents GENE simulations with a reduced velocity domain, $L_{v_\parallel }\times L_\mu = 4.5\times 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_\parallel } \times N_\mu = 6 \times 4$ grid does not yield unstable growth rates when the standard size of the velocity domain is considered.

3.2. 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 2(a) 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 figure 2(b), 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 $\eta _{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 2(b) also presents the convergence behaviour for two cases with larger numerical velocity dissipation, namely $\eta _{v}=0.01$ and $\eta _{v}=0.05$. In these cases, the moment approach converges more rapidly than with $\eta _{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_\parallel }\times N_\mu =8\times 4$ and $16\times 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 figure 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.

Figure 2. (a) Time traces of the heat flux for the CBC with comparison between the GM approach, GENE and the result of Dimits et al. (2000), $Q_x\approx 35$. Various velocity resolutions are used, in particular, different GM sets $(P,J)$ and numbers of points in the GENE velocity grid $N_{v_\parallel }\times N_\mu$. (b) Convergence in the saturated heat flux value with respect to the number of points representing the velocity space $N_p^v=(P+1)\times (J+1)$ for the GM approach and $N_p^v = N_{v_\parallel }\times N_\mu$ for GENE. The configuration space resolution is $N_x=128$, $N_y=64$ and $N_z=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 $L_v/2$ are obtained with a halved velocity domain size.

To further analyse the Hermite–Laguerre representation, we examine the amplitude of the individual GM, defined as $E_{pj}=\sum _{k_x,k_y,z}|N^{pj}(k_x,k_y,z)|$, time averaged during the quasi-steady saturated state shown in figure 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 detail, 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 for the differences observed between the GM energy spectra in the Hermite and Laguerre directions is twofold. First, the GM hierarchy (2.22) couples the Hermite modes mostly through the magnetic curvature and gradient drifts, connecting the $(p,j)$ GM with the $(p\pm 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\pm 1)$ GMs. This explains the oscillations observed in the Hermite modes ($p$ direction in figure 3). Second, we recall that the Laguerre polynomials are involved in the Bessel–Laguerre decomposition of (2.21), which converges non-monotonically with respect to the maximal Laguerre degree, explaining the non-monotonic convergence at the highest Laguerre modes in figure 3.

Figure 3. Amplitude of the GM, $E_{pj}=\sum _{k_x,k_y,z}|N^{pj}(k_x,k_y,z)|$, in the CBC, time averaged during the quasi-steady state.

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 figure 1a) corresponds to an increased saturated transport value (see figure 2b) and a high amplitude of the GM spectrum on figure 3. On the other hand, the $(8,4)$ basis presents a reduced amplitude for the high-degree GMs.

4. Collisionless Dimits shift

Building upon the study presented in § 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 § 3.

4.1. 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\leq k_y \leq 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 4(a) presents the results of a parameter scan involving the temperature gradient strength $\kappa _{T}$ and the number of GMs, having fixed $P=2J$, $\eta _{v}=0.001$ and all other parameters as in § 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)\gtrsim (12,6)$, the GM approach yields the threshold value from GK codes reported by Dimits et al. (2000), i.e. $\kappa _T\sim 4$. Furthermore, convergence is faster when the gradient level is larger, a trend consistent with Hoffmann et al. (2023). We recall here that this behaviour is a consequence of the relative increase of the importance of gradient drift terms (see (2.25) and (2.26)), compared with the parallel and perpendicular linear terms (see (2.23) and (2.24)), resulting in reduced coupling of high-degree moments.

Figure 4. (a) The ITG growth rate at $k_y\rho _s=0.3$ for different background temperature values $R_0/L_T$ and $(P,J)$, with $J=P/2$, $R_0/L_N=2.22$ and $\eta _{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 $R_0/L_T$ values between the ITG linear threshold and the non-zero $\chi$ values represents the Dimits shift.

4.2. 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_\parallel },N_\mu )=(8,4)$ and $(16,8)$ with parallel velocity and magnetic moment domain of size $L_{v_\parallel }\times L_\mu =4.5\times 1.5$. We also consider the $(N_{v_\parallel },N_\mu )=(32,16)$ and $(N_{v_\parallel },N_\mu )=(42,24)$ grids in the velocity domain for $L_{v_\parallel }\times L_\mu =9\times 3$. We measure the radial heat transport by evaluating the ion heat diffusivity, $\chi =\langle Q_x \rangle _t/(\kappa _T \kappa _N)$, with the heat flux averaged over time during the saturated phase of the nonlinear simulations.

Figure 4(b) 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 flow 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: $\eta _{v}=0.05$ (figure 5a), $\eta _{v}=0.01$ (figure 5b) and $\eta _{v}=0.001$ (figure 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 overestimate the transport values compared with Dimits et al. (2000) and GENE, for all dissipation levels and particularly for the lowest background gradient considered. This can be attributed to the overestimation of the linear growth rate, observed in particular at the lowest gradient levels. On the other hand, GENE linear results do not exhibit this overestimation at low resolution, which explains why the GENE results do not overestimate the transport at low resolution.

Figure 5. Heat diffusivity as a function of the temperature gradient strength $R/L_T$ obtained with the GM approach varying the intensity of the numerical velocity dissipation, $\eta _{v}=0.01$, $\eta _{v}=0.005$ and $\eta _{v}=0.001$ and obtained with GENE. The colours indicate the number of points in the velocity space, $N_p^v=(P+1)\times (J+1)$ with $(P,J)=(2,1)$, $(4,2)$, $(8,4)$, $(16,8)$ and $(30,15)$ for the GM approach and $N_p^v=N_{v\parallel }\times N_\mu =8\times 4$, $16\times 8$, $32\times 16$, $42\times 24$ for GENE. The results of GK and PIC simulations in Dimits et al. (2000) are reported in black.

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 $\eta _{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 compared 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.

5. 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, Sugama and Landau collision operators on the transport level. We first consider the same parameters as the CBC and set the temperature background level to $\kappa _T=5.3$. We vary the collision frequency parameter around $\nu \sim 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 $\kappa _T$, setting $\eta _{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 § 5), we first study the impact of our different collision models on the linear growth rates. These are presented in figure 6, with the collision frequency ranging from $\nu =0.05$ to $\nu =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 with 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 $\kappa _T=5.3$ for a set of collision frequencies, namely $0.005\leq \nu \leq 0.05$ with the $(P,J)=(16,8)$ polynomial basis. In figure 7(a), 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, Gürcan & Diamond 2015; Peeters et al. 2016; Hallenbert & Plunk 2022) as well as reduced turbulence modelling (Ivanov et al. 2020; Qi, Majda & Cerfon 2020; Ivanov, Schekochihin & Dorland 2022). Figure 7(b) 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, 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).

Figure 6. Collisional study of the ITG growth rate at $\kappa _T=5.3$ with the $(P,J)=(4,2)$ (blue), $(P,J)=(8,4)$ (orange) and $(P,J)=(16,8)$ (green) bases. Different collision operators are compared with collision frequencies $\nu =0.05$ (dashed) and $\nu =0.005$ (solid). The converged collisionless result are also shown (black circles).

Figure 7. Study of collision effects on the heat flux with the GM approach and the Dougherty (orange), Sugama (blue) and Landau (green) operators for $\kappa _T=5.3$ for $\nu =0.005$. The heat flux time traces (a) and the average heat diffusivity (b) are compared with the results from Lin et al. (1999), reported in grey. 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.

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 is 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 $\nu _{LDGK}=0.005$. As shown on figure 8, the collisionality tends to smooth out the transition between fully developed turbulence ($\kappa _T\gtrsim 6$) and zonal flows dominated saturated states ($\kappa _T\lesssim 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.

Figure 8. Heat diffusivity as a function of the temperature gradient strength $R/L_T$ obtained with the GM approach in the collisionless limit (blue circles) and for $\nu =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 $\kappa _T=6$ (dashed line).

6. Conclusions

In the present study, we report on a benchmark and convergence analysis of the GM approach in the local flux-tube $\delta 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\rho _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 with 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\times 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\times 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 the 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 ($\kappa _T=5.3$) reveals a minimal influence of collision on the ITG growth rate around a physical collision frequency of $\nu \sim 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 long-wavelength 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.