2.1. Nonlinear gyrokinetic model

Within the GK approach (Catto Reference Catto1978; Frieman & Chen Reference Frieman and Chen1982; Hazeltine & Meiss Reference Hazeltine and Meiss2003), 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 Reference Hazeltine and Meiss2003). 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 Reference Hazeltine and Meiss2003). 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. Reference Lapillonne, Brunner, Dannert, Jolliet, Marinoni, Villard, Görler, Jenko and Merz2009; Frei et al. Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2023)

(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 Reference Connor, Hastie and Taylor1978). 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. Reference Lapillonne, Brunner, Dannert, Jolliet, Marinoni, Villard, Görler, Jenko and Merz2009), we consider it here with the purpose of establishing direct comparisons with previous literature results (Lin et al. Reference Lin, Hahm, Lee, Tang and Diamond1999; Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000; Frei et al. Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2023).

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. Reference Giroud, Jachmich, Jacquet, Järvinen, Lerche, Rimini, Aho-Mantila, Aiba, Balboa and Belo2015). 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 Reference Grant and Feix1967; Madsen Reference Madsen2013; Manas et al. Reference Manas, Hornsby, Angioni, Camenen and Peeters2017; Adkins & Schekochihin Reference Adkins and Schekochihin2018; Staebler, Belli & Candy Reference Staebler, Belli and Candy2023)

(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 Reference Gradshteyn and Ryzhik2014).

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 Reference Ball, Brunner and Ajay2020). 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. Reference Frei, Jorge and Ricci2020). 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. Reference Neiser, Jenko, Carter, Schmitz, Told, Merlo, Banōn Navarro, Crandall, McKee and Yan2019; Hardman et al. Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022).

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. Reference Sugama, Watanabe and Nunami2009) and the GK Dougherty model (Dougherty Reference Dougherty1964). The details of these operators and their projection onto the Hermite–Laguerre basis can be found in Frei et al. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021). For an overview of the qualitative differences between the aforementioned collision operators we also refer to Hoffmann et al. (Reference Hoffmann, Frei and Ricci2023). 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.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. Reference Ball, Brunner and Ajay2020; Volčokas, Ball & Brunner Reference Volčokas, Ball and Brunner2023). 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. (Reference Beer, Cowley and Hammett1995) and Ball et al. (Reference Ball, Brunner and Ajay2020). 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. Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018) 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. (Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018).

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. (Reference Frei, Hoffmann and Ricci2022), Frei et al. (Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2023) and Hoffmann et al. (Reference Hoffmann, Frei and Ricci2023). 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 (Reference Zocco and Schekochihin2011) and Loureiro et al. (Reference Loureiro, Dorland, Fazendeiro, Kanekar, Mallet, Vilelas and Zocco2016) 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. (Reference Shukla, Hatch, Dorland and Michoski2022), which presents a method for formulating closures that learn from kinetic simulation data.

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. Reference Hoffmann, Frei and Ricci2023), 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) (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. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022; Frei et al. Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2023; Hoffmann et al. Reference Hoffmann, Frei and Ricci2023). 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. Reference Kotschenreuther, Rewoldt and Tang1995; Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008; Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011; Adkins, Ivanov & Schekochihin Reference Adkins, Ivanov and Schekochihin2023), 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). 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. Reference Pueschel, Dannert and Jenko2010), 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000), $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. Reference Hoffmann, Frei and Ricci2023). 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.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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000), i.e. $\kappa _T\sim 4$. Furthermore, convergence is faster when the gradient level is larger, a trend consistent with Hoffmann et al. (Reference Hoffmann, Frei and Ricci2023). 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. Reference Hoffmann, Frei and Ricci2023), 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) 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. Reference Hoffmann, Frei and Ricci2023). 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.