2.1. Gyrokinetic model in a Z-pinch configuration

We consider the GK approach (Catto Reference Catto1978; Frieman & Chen Reference Frieman and Chen1982; Hazeltine & Meiss Reference Hazeltine and Meiss2003) to study turbulence in a Z-pinch geometry. Using the standard $\delta f$ approach, we decompose $f_a$, the gyrocentre distribution function of species $a$ ($a=e$ for electrons and $a=i$ for ions), as the sum of a time-independent background Maxwellian component and a perturbation, $f_a = F_{aM} + \delta f_a$, where the Maxwellian distribution for a species $a$ is defined as $F_{aM} = N_a/({\rm \pi} ^{1/2}v_{\rm tha})^{3} \exp (-m_a v_\parallel ^2/2T_a-\mu B/T_a)$, with $\boldsymbol B=B \boldsymbol {b}$ the equilibrium magnetic field ($B = |\boldsymbol B|, \boldsymbol {b} = \boldsymbol B/B$), $N_a$ the equilibrium density, $T_a$ the equilibrium temperature, $m_a$ the particle mass, $\boldsymbol v=v_\parallel \boldsymbol {b} + \boldsymbol v_\perp$ the particle velocity, $\mu = m_a v_\perp ^2/B$ the particle magnetic moment and $v_{\rm tha}^2 = 2T_a/m_a$ the thermal velocity. We assume small fluctuations, $\delta f_a/F_{aM}\sim \varDelta \ll 1$, where the scaling parameter $\varDelta$ measures the perturbation amplitude relative to the background (Hazeltine & Meiss Reference Hazeltine and Meiss2003).

We focus here on the Fourier representation of the perturbed gyrocenter distribution function at the gyrocentre position $\boldsymbol R$ and a time $t$,

(2.1)\begin{equation} g_a(\boldsymbol k,v_\parallel,\mu,t):=\int \delta f_a(\boldsymbol R,v_\parallel,\mu,t) \, {\rm e}^{-{\rm i}\boldsymbol k \boldsymbol{\cdot} \boldsymbol R} \, {\mathrm d} \boldsymbol R, \end{equation}

using Fourier modes $\boldsymbol k = \boldsymbol k_\perp + k_\parallel \boldsymbol {b}$. The electrostatic GK Boltzmann equation determining the evolution of $g_a$ writes (Brizard & Hahm Reference Brizard and Hahm2007)

(2.2)\begin{equation} \partial_t g_a + \frac{1}{B}\{J_0\phi,g_a+F_{aM}\} + {\rm i}\omega_{Ba} h_a= \sum_{b}C_{ab}, \end{equation}

where we introduced the Poisson bracket operator, $\{\,f_1,f_2\}=\boldsymbol {b} \boldsymbol {\cdot } (\boldsymbol {\nabla } f_1 \times \boldsymbol {\nabla } f_2)$ for two generic fields $f_1,f_2$, to describe the effect of the background density and temperature gradients, and of the quadratic nonlinearities, of order $\varDelta ^2$, rising from the $\boldsymbol E\times \boldsymbol B$ drift. In (2.2) the magnetic drift frequency ${\rm i}\omega _{Ba}$ contains the magnetic curvature and gradient drifts, i.e.

(2.3)\begin{equation} \omega_{Ba} =\boldsymbol{b} \times \frac{1}{\varOmega_a}[{v_\parallel^2 (\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{b}+v_\perp^2 \boldsymbol{\nabla} B/B}]\boldsymbol{\cdot} \boldsymbol k_\perp,\end{equation}

where $\varOmega _a = q_a B/m_a$ is the cyclotron frequency with $q_a$ the particle charge. Following previous work (Ricci et al. Reference Ricci, Rogers and Dorland2006a; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), we assume $k_\parallel = 0$ in (2.2), and therefore we consider a two-dimensional domain that extends perpendicularly to the magnetic field line. The electrostatic potential $\phi$ is evaluated at the gyrocentre position through the gyroaveraging operator expressed, in Fourier space, with the zeroth-order Bessel function of the first kind, $J_0=J_0(b_a)$ with $b_a =|\boldsymbol k_\perp | |\boldsymbol v_\perp |/\varOmega _a$, containing FLR effects at all orders in $b_a$. Finally $C_{a,b}$ is the collision operator between species $a$ and $b$.

The electrostatic Poisson equation, in the quasi-neutrality limit, allows us to close the system by expressing the fluctuation of the electrostatic potential according to

(2.4)\begin{equation} \sum_a \frac{q_a^2}{N_a T_a}({1- \varGamma_0(b_{a,{\rm th}})})\phi = \sum_a q_a \int {\mathrm d} \boldsymbol v J_0 g_a, \end{equation}

where $b_{a,{\rm th}} = (k_\perp v_{\rm tha} /\varOmega _a)^2/2$ and $\varGamma _0(x) = I_0(x) \,{\rm e}^{-x}$ with $I_0$ the zeroth-order modified Bessel function of the first kind.

Equation (2.2) is now simplified considering the Z-pinch magnetic field and geometry. Using local field-aligned coordinates $(x,y,z)$, with $\boldsymbol {e}_x$ the radial, $\boldsymbol {e}_y$ the binormal and $\boldsymbol {e}_z$ the azimuthal directions, the Z-pinch magnetic field can be expressed as $\boldsymbol B = B \boldsymbol {e}_z$. The magnetic field presents a radial gradient, $\boldsymbol {\nabla } B/B = -1/L_B \boldsymbol {e}_x$, and curvature, $(\boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {b} = -1/L_B \boldsymbol {e}_x$, which are assumed constant within the flux-tube approach. The length $L_B$ denotes the distance between the flux tube and the symmetry axis of the Z-pinch (see figure 1). We also consider constant background density and temperature gradients, $\boldsymbol {\nabla } N_a/N_a = -1/L_N \boldsymbol {e}_x$ and $\boldsymbol {\nabla } T_a / T_a = -1/L_T\boldsymbol {e}_x$, for both electrons and ions. For comparison with common tokamak configuration, we note that, in the present geometry, all components of the metric tensor $g^{ij}$, for $i,j=x,y,z$, vanish except for $g^{xx}=g^{yy}=1$ and $g^{zz}=1/L_B^2$. In addition, one can express the Jacobian of the coordinate system as $J_{xyz}=L_B$ and the curvature operator as $[\boldsymbol {b}\times \boldsymbol {\nabla } B]\boldsymbol {\cdot }\boldsymbol {\nabla } = -B/L_B \partial _y$ where $\partial _y$ denotes the derivative in the $\boldsymbol {e}_y$ direction. We note that in the flux-tube framework, all background quantities ($B$,$N_a$ and $T_a$) and their associated gradients length ($L_B$, $L_N$ and $L_T$) are considered constant in time and in space.

Figure 1. Illustration of the Z-pinch magnetic geometry considered and the simulated perpendicular plane (grey area). The field-aligned coordinate system and a magnetic field line $\boldsymbol B$ (blue arrow) are depicted. We also indicate the direction of the density and temperature equilibrium gradients, $\boldsymbol {\nabla } N$ and $\boldsymbol {\nabla } T$, in addition to the magnetic equilibrium gradient and curvature, $\boldsymbol {\nabla } B$ and $\boldsymbol b\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol b$, respectively (orange arrow). The symmetry axis of the cylinder is represented by the dashed-dotted line and $L_B$ denotes the distance between the cylinder axis of symmetry and the flux tube.

Throughout the rest of this work, we use the following dimensionless units. The dimensionless parallel and perpendicular velocity coordinates are defined by $s_{\parallel a} = v_\parallel /v_{\rm tha}$ and $x_a=\mu B /T_a$, respectively. The perpendicular spatial scales are normalized to the sound Larmor radius $\rho _{s}=c_s/\varOmega _i$, with $c_s=\sqrt {T_e/m_i}$ the sound speed. Time is normalized to $L_B/c_s$. The electrostatic potential is normalized to $T_e/e$ with $e$ the elementary charge, which allows us to define the normalized particle charge, $z_a=q_a/e$, as well. We define the temperature and mass ratio $\tau _a=T_a/T_e$ and $\sigma _a=\sqrt {m_a/m_i}$, respectively, and we introduce the dimensionless density gradient drive, $\kappa _N=L_B/L_N$, the dimensionless temperature gradient drive, $\kappa _T=L_B/L_T$, and their ratio, $\eta = |\boldsymbol {\nabla } \ln T|/|\boldsymbol {\nabla } \ln N|$. It is worth noting that the flux-tube limit in a Z-pinch is valid for $\rho _s/L_B\ll 1$.

Considering purely perpendicular Fourier modes $\boldsymbol k = \boldsymbol k_\perp = k_x \boldsymbol {e}_x + k_y \boldsymbol {e}_y$, the GK equation for $g_a$, (2.2), writes

(2.5)\begin{align} \partial_t g_a +\{J_0\phi,g_a\} - \frac{\tau_a}{z_a}[{2s_{{\parallel} a}^2+x_a}] \, {\rm i} k_y h_a +\left[\kappa_N +\kappa_T\left(s_{{\parallel} a}^2+x_a-\frac{3}{2}\right)\right]\,{\rm i}k_y J_0 \phi = \sum_{b}C_{ab}, \end{align}

where we introduced the non-adiabatic part of the distribution perturbed function

(2.6)\begin{equation} h_a(\boldsymbol k, s_{{\parallel} a}, x_a, t) = g_a(\boldsymbol k, s_{{\parallel} a}, x_a, t) + z_a/\tau_a J_0 \phi(\boldsymbol k, t). \end{equation}

The Poisson bracket in the Z-pinch geometry writes as $\{f_1,f_2\}=\partial _x f_1 \partial _y f_2 - \partial _y f_1 \partial _x f_2$ in real space. In the Fourier, this yields a convolution expressed as

(2.7)\begin{equation} \{f_1,f_2\}=\sum_{k_x',k_y'} k_x (k_y-k_y') f_1[\boldsymbol k-\boldsymbol k']f_2[\boldsymbol k] - k_y (k_x-k_x') f_1[\boldsymbol k-\boldsymbol k'] f_2[\boldsymbol k]. \end{equation}

Finally, we close our system with the dimensionless Poisson equation, i.e.

(2.8)\begin{equation} \sum_a \frac{z_a^2}{\tau_a}({1- \varGamma_0(b_{a,{\rm th}})})\phi = \sum_a z_a \int {\mathrm d} s_{{\parallel} a} \,{\mathrm d} x_a J_0 g_a.\end{equation}

2.2. Nonlinear gyromoment hierarchy

In order to solve (2.5) by using the gyromoment framework, we expand the distribution function on a Hermite–Laguerre polynomial basis (Jorge, Ricci & Loureiro Reference Jorge, Ricci and Loureiro2017; Frei et al. Reference Frei, Jorge and Ricci2020), i.e.

(2.9)\begin{equation} g_a(\boldsymbol k, s_{{\parallel} a}, x_a, t) = \sum_{p,j}N_a^{pj}(\boldsymbol k,t)H_p(s_{{\parallel} a})L_j(x_a)F_{aM}(s_{{\parallel} a},x_a). \end{equation}

In (2.9), we introduce the gyromoment of order $(\,p,j)$, i.e. the basis coefficient

(2.10)\begin{equation} N_a^{pj}(\boldsymbol k,t)=\int_0^\infty {\mathrm d} x_a \int_{-\infty}^\infty {\mathrm d} s_{{\parallel} a} g_a(\boldsymbol k, s_{{\parallel} a}, x_a, t) H_p(s_{{\parallel} a}) L_j(x_a), \end{equation}

where

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

and

(2.12)\begin{equation} L_j(x_a)=\frac{{\rm e}^{x_a}}{j!}\frac{{\rm d}^j}{{\rm d}x_a^j}x_a^j \, {\rm e}^{-x_a},\end{equation}

are the physicist's Hermite polynomial of order $p$ and the Laguerre polynomial of order $j$, respectively (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014). The Hermite polynomials of (2.11) are normalized such that $\int _{-\infty }^\infty {\mathrm d} s_{\parallel a} H_p H_{p'} \, {\rm e}^{-s_{\parallel a} ^2}=\delta _{pp'}$ where $\delta _{pp'}$ denotes the Kronecker delta. Similarly, the Laguerre polynomials satisfy the orthogonality relation ${\int _{0}^\infty {\mathrm d} x_a L_j L_{j'} \, {\rm e}^{-x_a}=\delta _{jj'}}$. The use of Hermite polynomial projection is common in the literature, particularly for projecting the one-dimensional velocity space Vlasov–Poisson system (Armstrong Reference Armstrong1967; Grant & Feix Reference Grant and Feix1967; Joyce, Knorr & Meier Reference Joyce, Knorr and Meier1971; Gibelli & Shizgal Reference Gibelli and Shizgal2006; Parker & Dellar Reference Parker and Dellar2015). On the other hand, Laguerre polynomials are not as frequently used in plasma physics compared with Hermite polynomials. Aside from spanning fluid equations (Manas et al. Reference Manas, Hornsby, Angioni, Camenen and Peeters2017), their main application is in expressing collision models in a spectral framework (Brunner, Valeo & Krommes Reference Brunner, Valeo and Krommes2000; Belli & Candy Reference Belli and Candy2012).

We now project the Boltzmann GK equation, (2.5), onto the Hermite–Laguerre basis. We expand the Bessel function of the first kind in terms of Laguerre polynomials as

(2.13)\begin{equation} J_0 = J_0(\sqrt{l_a x_a}) = \sum_{n=0}^\infty \mathcal{K}_n(l_a)L_n(x_a),\end{equation}

with the kernel functions $\mathcal {K}_n(l_a)=l_a^n \, {\rm e}^{-l_a}/n!$, being $l_a = \sigma _a^2\tau _ak_\perp ^2/2$ (Frei et al. Reference Frei, Jorge and Ricci2020). The projection of (2.5) yields the gyromoment nonlinear hierarchy in a Z-pinch configuration, which can be expressed as

(2.14)\begin{equation} \partial_t N_a^{pj} + \mathcal{S}_a^{pj} + \mathcal{M}_a^{pj} + \mathcal{D}_a^{pj} = \mathcal{C}_a^{pj}, \end{equation}

where the term related to the magnetic gradient and curvature drifts yields

(2.15)\begin{align} \mathcal{M}_a^{pj} & =-\frac{\tau_a}{z_a} \, {\rm i} k_y [{\sqrt{(\,p+1)(\,p+2)} n_a^{p+2,j} + (2p+1)n_a^{p,j} + \sqrt{p(\,p-1)}n_a^{p-2,j}}]\nonumber\\ & \quad - \frac{\tau_a}{z_a} \, {\rm i} k_y [{(2j+1)n_a^{pj} - (\,j+1)n_a^{p,j+1}-jn_a^{p,j-1}}]. \end{align}

The term related to the density and temperature gradients writes

(2.16)\begin{equation} \mathcal{D}_a^{pj}= \kappa_N \,{\rm i}k_y\phi[\mathcal{K}_j\delta_{p0}+\eta\mathcal{K}_j\frac{\sqrt{2}}{2}\delta_{p2} + \eta(2j\mathcal{K}_j-[\,j+1]\mathcal{K}_{j+1}-j\mathcal{K}_{j-1})\delta_{p0}].\end{equation}

In (2.15) and (2.16), we introduce the non-adiabatic gyromoments $n_a^{pj}(\boldsymbol k,t)=N_a^{pj}+z_a/\tau _a \mathcal {K}_j \phi \delta _{p0}$.

The Hermite polynomial product rule, $s_{\parallel a} H_p = \sqrt {(\,p+1)/2}H_{p+1} + \sqrt {p/2}H_{p-1}$, and the Laguerre polynomial product rule, $x_a L_j = (2j+1)L_j-(\,j+1)L_{j+1}-jL_{j-1},$ are used to deduce (2.15) and (2.16).

The nonlinear term related to the $\boldsymbol E\times \boldsymbol B$ drift is expressed in terms of gyromoments by using the Bessel–Laguerre decomposition, (2.13), and the Poisson bracket, (2.7), which yields

(2.17)\begin{equation} \mathcal{S}_a^{pj} = \sum_{n=0}^{\infty}\left\{\mathcal{K}_n\phi,\sum_{s=0}^{n+j}d_{njs} N_a^{ps}\right\}. \end{equation}

To obtain (2.17), we expressed the product of two Laguerre polynomials as a sum of single polynomials using the identity

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

with

(2.19)\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!} \left(\begin{array}{@{}c@{}}{n}\\{n_1}\end{array}\right)\left(\begin{array}{@{}c@{}}{j}\\{j_1}\end{array}\right)\left(\begin{array}{@{}c@{}}{s}\\{s_1} \end{array}\right).\end{equation}

This choice differs from the representation in the GX code that evaluates the Laguerre product with a pseudo-spectral algorithm in the velocity space (Mandell et al. Reference Mandell, Dorland and Landreman2018a, Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022).

The Poisson equation, (2.8), is also projected onto the Hermite–Laguerre basis. This yields (Frei et al. Reference Frei, Jorge and Ricci2020)

(2.20)\begin{equation} \left[\sum_a \frac{z_a^2}{\tau_a}\left(1-\sum_{n=0}^{\infty}\mathcal{K}^2_n\right)\right]\phi = \sum_a z_a \sum_{n=0}^{\infty}\mathcal{K}_n N_a^{0n},\end{equation}

where the quasi-neutrality approximation is used, i.e. $(k_\perp \lambda _D)^2\ll 1$ with $\lambda _D$ the Debye length. In the collisionless limit ($\mathcal {C}_a^{pj}=0$), the gyromoment hierarchy, (2.14), combined with the Poisson equation, (2.20), can be considered as an extension of the gyrofluid model to an arbitrary number of moments. In Appendix A, we demonstrate that the formerly derived gyrofluid model in Brizard (Reference Brizard1992) can be retrieved by properly truncating the collisionless gyromoment hierarchy. We also show how an extended Hasegawa–Mima model (Hasegawa & Mima Reference Hasegawa and Mima1978; Dewhurst, Hnat & Dendy Reference Dewhurst, Hnat and Dendy2009) can be obtained. It is worth noting that the gyromoment approach presented here conserves free energy in the collisionless limit and when driving gradients are neglected, even when closure by truncation is used (Mandell et al. Reference Mandell, Dorland and Landreman2018a).

Finally, we note that we characterize the turbulent transport in a Z-pinch by considering the dimensionless ion particle flux, $\boldsymbol \varGamma = n_i \boldsymbol v_{E\times B}$, with $\boldsymbol v_{E\times B}=-\boldsymbol {\nabla } \phi \times \boldsymbol {b}$ the $\boldsymbol E\times \boldsymbol B$ velocity and $n_i=\sum _{n=0}^\infty \mathcal {K}_n N_i^{0n}$ the ion particle density perturbation. In the following, we analyse the time series of the spatially averaged radial ion particle flux, $\varGamma _x(t) = \langle \boldsymbol \varGamma \boldsymbol {\cdot } \boldsymbol {e}_x \rangle _{xy}$, which can be expressed, using the Fourier modes of the gyromoments, as

(2.21)\begin{equation} \varGamma_x(t) =\sum_{k_x,k_y} ({\rm i}k_y\phi)^* \sum_{n=0}^{\infty}\mathcal{K}_n N_i^{0n}. \end{equation}

The saturated radial particle transport, $\varGamma _x^{\infty }$, is analysed by evaluating the convergence of the quantity $\overline {\varGamma _x}(t) = \int _{t_0}^{t}\varGamma _x(t')\, {\rm d} t'/(t-t_0)$ as $t$ increases, considering $t_0$ sufficiently large that the initial transient present in the simulation is not considered. The value of $\bar \varGamma _x(t)$ provides an estimate for the saturated transport level, $\varGamma _x^{\infty } = {\lim }_{t\rightarrow \infty }\bar \varGamma _x(t)$.

2.3. Linear collision operators

The $\mathcal {C}_a^{pj}$ term in (2.14) represents the effect of collisions through the projection of a collision operator model onto the Hermite–Laguerre basis. Any linearized Fokker–Planck collision operator can be written as the sum of a test part $C^T_{ab}$ and a field part $C^F_{ab}$, i.e. $C_{ab} = C^T_{ab} +C^F_{ab}$ (Helander & Sigmar Reference Helander and Sigmar2002; Hazeltine & Meiss Reference Hazeltine and Meiss2003), with $C^T_{ab} = C(\,f_a,F_{bM})$ and $C^F_{ab} = C(F_{aM},f_b)$ for any species $a$ and $b$.

We consider here the Coulomb, Sugama, Lorentz and Dougherty operators. For the case of the Coulomb collision operator, we introduce the Rosenbluth potentials, $H(\,f)=2\int {\rm d}^3 v' f(\boldsymbol v')/|\boldsymbol v-\boldsymbol v'|$ and $G(\,f)=\int {\rm d}^3 v' |\boldsymbol v - \boldsymbol v'| f(\boldsymbol v')$, as well as the phase space coordinates $(\boldsymbol r, v, \xi,\theta )$, where $\boldsymbol r$ denotes the particle position, $v$ the magnitude of its velocity, $\xi =v_\parallel /v$ the pitch angle and $\theta$ the gyro-angle. In this framework, the test part of the Coulomb collision operator can be expressed as (Frei et al. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021)

(2.22)\begin{align} C_{ab}^T & = \frac{m_a\nu_{ab}}{N_b}\left\{ \left[\frac{2}{v^2}G(F_{bM} ) + \left(1 - \frac{m_a}{m_b}\partial_v H(F_{bM})\right)\right] \partial_v f_a \right.\nonumber\\ & \quad \left. - \frac{1}{v^3} \partial_v G(F_{bM}) \mathcal{L}^2 f_a + \partial^2_v G(F_{bM})\partial_v^2 f_a + \frac{m_a}{m_b}8{\rm \pi} F_{bM} f_a \right\}, \end{align}

where $\mathcal {L}f=\partial _\xi [(1-\xi ^2)\partial _\xi f] + \partial ^2_\theta f/(1-\xi ^2)$ is the pitch-angle operator. On the other hand, the field part yields

(2.23)\begin{align} C_{ab}^F=\frac{2\nu_{ab}v_{\rm tha} F_{aM}}{N_b} \left[2\frac{v^2}{v_{\rm tha}^2}\partial_v^2 G(\,f_b) - H(\,f_b) - \left(1-\frac{m_a}{m_b}\right)v\partial_v H(\,f_b) + \frac{m_a}{m_b}4{\rm \pi} v_{\rm tha}^2 f_b\right].\end{align}

It is worth noting that the computation of the field term is particularly costly because of the velocity integrals of the perturbed distribution function contained in the Rosenbluth potentials. However, the projection of this operator on the Hermite–Laguerre basis enables the expression of these integrals as a linear combination of gyromoments. In this work, the Coulomb operator refers to the gyro-averaged version of the linearized Fokker–Planck collision operator, (2.22) and (2.23).

The Sugama collision model (Sugama et al. Reference Sugama, Watanabe and Nunami2009) is a multi-species generalization of the Abel operator (Abel et al. Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008). While the Sugama operator considers the test part of the Fokker–Planck operator, (2.22), which includes pitch-angle scattering and energy diffusion, the field term is replaced by an ad hoc term derived from a fluid approach to conserve particle, momentum, energy and satisfy the H-theorem.

The Lorentz model considers like-particle collisions, i.e. $a=b$, and it is based on the small mass ratio limit, which simplifies the test part of (2.22) to the pitch-angle operator term only. The field part is adapted to conserve particles, momentum and energy. Ricci et al. (Reference Ricci, Rogers and Dorland2006a) observe that this operator does not provide sufficient damping to avoid the use of artificial dissipation in nonlinear simulations. This is in contrast to the Sugama and Abel operators.

Finally, the Dougherty model consists of kinetic and spatial second-order diffusion terms (Lenard & Bernstein Reference Lenard and Bernstein1958) with corrections involving the density, velocity and temperature fluid moments in order to conserve particle, momentum and energy. The details of the Dougherty, Sugama, Lorentz and Coulomb GK operators as well as their projection onto the Hermite–Laguerre basis can be found in Frei et al. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021). We set the intensity of the collisions through the normalized ion–ion collision frequency $\nu$. The collision frequencies among the different species are thus given by $\nu _{ii}=\nu$, $\nu _{ee}=\sigma _e \tau _e^{3/2}\nu$, $\nu _{ei}=\nu$ and $\nu _{ie}=\sigma _e \tau _e^{3/2}\nu$.

2.4. Numerical approach

To solve (2.14) numerically, we evolve a finite set of gyromoments $N_a^{pj}(\boldsymbol k,t)$ with $0\leq p \leq P$ and $0\leq j \leq J$ and consider the Fourier modes with $k_x=m{\rm \Delta} k_x$, with ${0\leq m\leq M}$, and $k_y=n{\rm \Delta} k_y$, with $-N/2+1\leq n \leq N/2$, using a standard explicit fourth-order Runge–Kutta time-stepping scheme. In the Z-pinch geometry, the gyromoments hierarchy decouples odd and even Hermite gyromoments, which is a consequence of the $k_\parallel =0$ assumption. This allows us to evolve only the even gyromoments $N_a^{pj}$, with $p= {2l,l\in \mathbb {N}}$. The hierarchy is closed by using a simple truncation, i.e. $N_a^{p,j}=0$ for all $p>P$ or $j>J$. In the following, we denote this truncated gyromoment set as a $(P,J)$ basis. The use and analysis of more advanced closure schemes, e.g. 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), are left for future work.

Focusing on the nonlinear term in (2.17), we first observe that any truncation of the sum over $s$ must be avoided in order to prevent polynomial aliasing. Hence, to guarantee the exact Laguerre product identity in (2.18), we truncate the sum over $n$ in (2.17) to $n\leq J-j$. Second, we note that the computation of the $d_{njs}$ coefficients is challenging since they involve sums and differences of large numbers. To avoid the overflow of the floating point representation, we use an arbitrary precision library for our calculations (Smith Reference Smith1991). Finally, we note that the convolutions in Fourier space are treated with a conventional pseudo-spectral method, i.e. the backward fast Fourier transform (Frigo & Johnson Reference Frigo and Johnson2005) of the fields to convolve, the multiplication in real space and the forward fast Fourier transform of the result, including the usual 2/3 Orszag rule for anti-aliasing (Orszag Reference Orszag1971).

Regarding the collision operators, the Dougherty operator, which has a light computational cost, is directly implemented in the gyromoment hierarchy. On the other hand, the evaluation of the Sugama, Lorentz pitch-angle and Coulomb collision terms is reduced to a four-dimensional matrix-vector operation, i.e. the $pj$th collision term is written as $C_{a}^{pj}=\sum _b\sum _{p'=0}^{P_b}\sum _{j'=0}^{J_b}\mathcal {C}_{ab}^{pj,p'j'} N_{b}^{p'j'}$ with a precomputed collision matrix $\mathcal {C}_{ab}^{pj,p'j'}$ of size $(P_a\times J_a)\times (P_{b}\times J_{b})$. The projection of the collision operators on the Hermite–Laguerre basis and the details of the computation of the matrix coefficients for each collision operator considered in the present work can be found in Frei et al. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021). It is worth noting that the GK corrections create a $k_\perp$ dependence of the matrix coefficients, i.e. $\mathcal {C}_{ab}^{pj,p'j'}=\mathcal {C}_{ab}^{pj,p'j'}(k_\perp )$, that calls for the precomputation of the coefficients for each $k_\perp$ present in the simulations. The computational cost of evaluating the GK matrix coefficients increases with $k_\perp$. Indeed, at large $k_\perp$, accurate FLR effects ask for larger bounds in the truncated sums used to approximate Bessel functions and basis transformations (see Frei et al. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021) for more details). We ensured the convergence of our matrix evaluation by analysing the eigenvalue spectrum and the matrix symmetry.

3.1. Entropy mode instability

At density gradients below the magneto-hydrodynamic (MHD) interchange instability threshold, a small-scale non-MHD instability, the entropy mode, can be destabilized in the Z-pinch configuration. The region of stability of the entropy mode is presented in Ricci et al. (Reference Ricci, Rogers, Dorland and Barnes2006b). Our analysis focuses on density gradient values $1.6\leq \kappa _N \leq 2.5$, while the temperature density gradient ratio is constant, $\eta =0.25$. This parameter encompasses an unstable region of the entropy mode, which extends, indeed, from $\kappa _N\simeq 2.5$, where the ideal MHD interchange mode is destabilized, to $\kappa _N\simeq 1.6$, which is close the analytical stability limit found by Ricci et al. (Reference Ricci, Rogers, Dorland and Barnes2006b), that is $\kappa _N={\rm \pi} /2$ for $\eta = 0$.

We start the analysis of the collisionless case by focusing on the linear growth rate of the entropy mode. The entropy mode instability growth rate is obtained by solving the initial value problem associated with the gyromoment hierarchy, (2.14), coupled to the Poisson equation, (2.20), and where the nonlinear terms, developed in (2.17), are neglected. In particular, we evolve $\phi$ and $N_a^{pj}$ as a function of $k_y$ modes, setting $k_x=0$ where the entropy mode growth rate peaks (Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006b). We compute the growth rates, $\gamma (k_y)$, by fitting the slope of the time evolution of $\ln |\phi _{k_y}|$ over a time window. The convergence is tested by checking that the results are independent of the size of the time window.

Considering the results presented on figure 2, we first note that the gyromoment approach retrieves the converged results obtained with GENE, given a sufficiently large polynomial basis. The convergence properties of the gyromoment model depend on the strength of the gradients and improve at steep gradients, confirming previous results obtained for the slab (Frei et al. Reference Frei, Hoffmann and Ricci2022a) and the toroidal ITG instability (Frei et al. Reference Frei, Hoffmann, Ricci, Brunner and Tecchiolli2022b). Second, the results obtained with a number of polynomials below convergence show a stabilization of the high $k_y$ tail of the entropy mode and a larger peak growth rate. Third, it is worth noting that, independently of the polynomial resolution, the growth rates obtained with a small number of polynomials agree with the converged results in the long-wavelength limit, $k_y\ll 1$, highlighting the fact that the gyromoment method retrieves the fluid limit, even when a small set of gyromoments is used.

Figure 2. Growth rates of the linear entropy mode in the collisionless case ($\nu =0$) for three different drive values, $\kappa _N=1.6$ (a), $\kappa _N=2.0$ (b) and $\kappa _N=2.5$ (c), keeping $\eta =0.25$. The growth rates are obtained with GENE with $N_{v_\parallel }=32$ and $N_\mu = 16$ velocity grid points (black diamonds) and different gyromoment sets: $(4,2)$ (blue), $(10,5)$ (red), $(20,10)$ (yellow) and $(30,15)$ (purple).

3.2. Nonlinear collisionless simulations

Let us now consider the nonlinear case by including the $E\times B$ term in (2.17). GENE results are used to benchmark our implementation (the GENE simulations presented here closely recall those by Hallenbert & Plunk Reference Hallenbert and Plunk2022). We focus on three values of the density background gradient, $\kappa _N = 1.6$, $2.0$ and $2.5$, with $\eta = 0.25$ and $\nu =0$. The system is evolved in a periodic box of dimensions $L_x\times L_y=120\times 80$, for the lowest gradient value, and $L_x\times L_y=400\times 240$, for the highest gradient value. In terms of spatial resolution, we consider a Fourier grid with $N=128$ and $M=32$ Fourier modes along the $x$ and $y$ directions, respectively, except for the steepest gradient case where we increase the resolution to $N=256$ and $M=128$ in order to reduce the need of artificial numerical dissipation. The velocity space is represented by the Hermite–Laguerre basis $(P,J)=(4,2)$ extended up to $(P,J)=(20,10)$ at the lowest gradient. GENE results are obtained using the same Fourier modes as the gyromoment simulation and a velocity grid resolution of $N_{v_\parallel }\times N_\mu =32\times 12$ points for the $(v_\parallel,\mu )$ velocity space in a box of dimension $L_{v_\parallel } \times L_\mu = 6\times 4$. This ensures convergence of GENE results and that the results obtained by Hallenbert & Plunk (Reference Hallenbert and Plunk2022) are retrieved. It is worth noting that Hallenbert & Plunk (Reference Hallenbert and Plunk2022) present the velocity space resolution we use as the minimum necessary not to compromise key results. We confirm their claim by observing spurious predator–prey cycles when running lower resolution simulations at $\kappa _N=1.6$.

When running the collisionless cases, GENE uses a kinetic artificial diffusion term, $\nu _v ({\rm \Delta} v_\parallel /2)^4\partial _{v\parallel }^4 g_a$, with the diffusion parameter fixed to $\nu _v=0.2$ (Pueschel, Dannert & Jenko Reference Pueschel, Dannert and Jenko2010). Both codes use a spatial fourth-order hyperdiffusion term in both perpendicular directions $\mu _{\rm HD} (k/k_{\mathrm {max}})^4$, with $0.5\leq \mu _{\rm HD}\leq 5.0$, adjusted on the drive level in order to avoid energy pile up without compromising the accuracy of results. For the intermediate values of the equilibrium gradient strength, we perform two simulations with GENE. The first simulation considers a constant level of numerical diffusion, while the second takes advantage of the adaptive numerical diffusion feature in GENE, as described in Hallenbert & Plunk (Reference Hallenbert and Plunk2022). This is done to ensure that the effect of the adaptive diffusion feature is not significant. This test allows us to confirm that the level of transport is resilient to spatial hyperdiffusion.

While a comparison of the computational cost of the two approaches is not straightforward, we note that the number of gyromoments evolved is given by $N_{P,J}=(P/2+1)\times (J+1)$ (we take into account that only the even $p$ gyromoments are evolved in the Z-pinch geometry). Therefore, 9 and 25 gyromoments are evolved in the $(P,J)=(4,2)$ and $(20,10)$ simulations, respectively. This compares with the, approximately, $10^2$ velocity grid points used by GENE.

We now focus on the quasi-steady turbulent state that is established after an initial transient following the initialization of the simulation. In particular, we measure the saturated time-averaged turbulent transport level, $\varGamma _x^{\infty }$. We observe that the gyromoment approach retrieves the saturated turbulent transport level obtained by GENE for all gradient values, given a sufficient number of gyromoments, over four orders of magnitudes (see figure 3). As for the linear case, faster convergence with the number of gyromoments is observed in the case of the strongest gradient, with a set of $(P,J)=(4,2)$ gyromoments being sufficient for convergence. This result might be surprising, considering the linear growth rate obtained with the same gyromoment resolution, significantly broader and showing a higher peak value than the converged value (see figure 2). Even the results obtained with $(P,J)=(4,2)$ at the lowest gradient are surprisingly accurate when considering the accuracy of the linear growth rate.

Figure 3. Comparison of the time-averaged transport level $\varGamma _x^{\infty }=\langle \varGamma _x\rangle _t$ obtained with the gyromoment (GM) approach for $(P,J)=(4,2)$ (blue), $(10,5)$ (red) and $(20,10)$ (yellow) and GENE, $\eta =0.25$. The time traces are presented on figure 5.

One can explain the faster convergence of the nonlinear simulations with respect to the evolution of the linear growth rate by considering the Fourier spectrum of the radial particle transport, $\langle \boldsymbol |\boldsymbol \varGamma \boldsymbol {\cdot } \boldsymbol {e}_x|\rangle _t$, at $k_x=0$ (see figure 4). It is found that transport is driven by fluctuations that occur on scale lengths that are larger than the ones at the peak growth rate of the entropy mode (see figure 2). The peak of the transport spectrum shifts towards smaller wavelengths when the density gradient is reduced, in good agreement with the transport scaling, $\varGamma _x\sim \gamma ^2/k^3$, derived in § 4. In addition, the turbulent quasi-steady state at low driving gradients is dominated by ZF which result from the growth of a KHI rising from $E\times B$ shear flow produced by the primary instability. The growth rate of the KHI typically peaks at wavelengths that are twice as long as the primary instability (Rogers & Dorland Reference Rogers and Dorland2005), thus pushing the dynamics towards larger spatial scales, where convergence of the gyromoment approach is achieved with a smaller number of gyromoments.

Figure 4. Spectrum of the radial particle flux, $\langle |\varGamma _x(k_x=0,k_y)|\rangle _t$ for the highest resolution simulations presented in figure 6 for $\kappa _N=2.5$ (solid squares), $\kappa _N=2.0$ (dashed diamonds) and $\kappa _N=1.6$ (dotted circles), $\eta =0.25$.

Comparing the convergence of the entropy mode growth rate (figure 2), the convergence of saturated transport level (figure 3) and the spectrum of the radial particle transport (figure 4), one can infer that the gyromoment simulations yield an accurate nonlinear transport level when the linear growth rate of the entropy mode is converged at the wavenumber of the transport spectrum peak. For example, considering the $\kappa _N=2.0$ case, we note that the transport spectrum peaks at $k_y\simeq 0.3$. The $(P,J)=(4,2)$ gyromoment result provides accurate growth rates for $k_y\lesssim 0.2$, thus yielding an inaccurate transport level. On the other hand, the $(P,J)=(10,5)$ gyromoment set is linearly accurate for $k_y\lesssim 0.5$, which explains the correct saturated transport result.

For a finer analysis of our simulation results, we study the time dependence of the turbulent transport for the three equilibrium gradient values (see figure 5). For instance, we note a negligible variation of the transport level with respect to the hyperdiffusion parameter, which mostly affects small-scale fluctuations, confirming Ricci et al. (Reference Ricci, Rogers and Dorland2006a) and Hallenbert & Plunk (Reference Hallenbert and Plunk2022). At a large gradient level, $\kappa _N=2.5$, the gyromoment approach qualitatively and quantitatively agrees with GENE, showing an approximately constant transport. The analysis of the turbulent eddies shows fully developed turbulence, with a negligible role of ZF (see figure 6). At the intermediate gradient value, $\kappa _N=2.0$, time intervals characterized by a high turbulent transport level ($\varGamma _x\sim 1$) alternate with quiescent periods ($\varGamma _x\ll 1$), as shown in figure 6. The $(P,J)=(10,5)$ simulation is in good agreement with GENE results, while the $(P,J)=(4,2)$ results underestimate the average $\varGamma _x$ value because of longer low-transport intervals and lower burst level. However, our numerical tests show that the level of agreement between the gyromoment approach and GENE is within an uncertainty similar to the one related to the use of a constant hyperdiffusion or an adaptive numerical diffusion algorithm in GENE. Finally, at the lowest gradient value considered, $\kappa _N=1.6$, the system is dominated by strong ZF that quench the turbulence reducing drastically the transport (see figure 6). As expected, the gyromoment method shows the largest discrepancies with respect to GENE in this case. GENE simulation results in transport with small amplitude fluctuations occurring on long time scales around a plateau value, $\varGamma _x^{\infty }\sim 10^{-2}$. On the other hand, the gyromoment approach shows bursts related to the damping of the ZF, for both the $(P,J)=(4,2)$ and $(10,5)$ resolutions. Hence, even though the $(P,J)=(10,5)$ simulation results in an averaged transport level similar to GENE, an accurate description of the turbulent dynamics requires a larger number of gyromoments. This is demonstrated by a $(P,J)=(20,10)$ simulation (see yellow line in figure 5), which agrees better with GENE results and does not produce the spurious bursts observed when a lower number of gyromoments is used.

Figure 5. Radial particle transport $\varGamma _x(t)$ (see (2.21)) from our nonlinear simulations. GENE results are obtained with constant $\mu _{\rm HD}$ (black) and $\mu _{\rm HD}$ set by an adaptive hyperdiffusion algorithm (grey). The gyromoment results are shown for $(P,J)=(4,2)$ (blue), $(P,J)=(10,5)$ (red) and $(P,J)=(20,10)$ (yellow). In all cases $\eta =0.25$.

Figure 6. Snapshots of the electrostatic potential (a,c,e) and the charge density $n_i-n_e$ (b,d,f) in the collisionless case at the three drive values considered, i.e. $\kappa _N=2.5$ (a,b) with $(P,J) = (4,2)$, $\kappa _N=2.0$ with $(P,J) = (10,5)$ (c,d) and $\kappa _N = 1.6$ with $(P,J) = (20,10)$ (e,f). In all cases, $\eta =0.25$.

We note that bursts can also be obtained with GENE by reducing the $(v_\parallel,\mu )$ velocity grid resolution to $16\times 8$, keeping $\nu _v=0.2$. Bursts are also obtained with a $32\times 16$ resolution when the velocity diffusion parameter is increased by a factor of $16$, i.e. $\nu _v = 3.2$, which ensures the same level of dissipation as in the coarser velocity resolution case. Thus, predator–prey cycles appear when a large level of diffusion is present in the velocity space. This diffusion can also be introduced through simple collision models, such as the Lenard–Bernstein operator (Lenard & Bernstein Reference Lenard and Bernstein1958). Our results thus demonstrate that the effect of using a reduced number of gyromoments is comparable to the presence of diffusion in velocity space, with the level of diffusion that depends on the highest gyromoment considered.

Since the representation of the velocity dependence of the distribution functions differs fundamentally between gyromoments and continuum approaches, we compare the time-averaged velocity distribution functions obtained by the gyromoments and GENE codes. Within the gyromoments method, one can reconstruct the distribution function by using the gyromoments as coefficients of the Hermite–Laguerre basis. This yields the averaged velocity distribution

(3.1)\begin{equation} g_{v,a}(s_{{\parallel} a},x_a,t) = \sum_{p=0}^P \sum_{j=0}^J \langle N_a^{pj}(\boldsymbol k, t) \rangle_{k_x,k_y} H_p(s_{{\parallel} a}) L_j(x_a) F_{aM}. \end{equation}

The results are presented in figures 7 for the ion distribution functions, considering $\kappa _N=1.6$ and $2.5$. As for the transport properties, the agreement of the distribution functions between both codes depends on the gradient value. At all gradient values considered, the $(P,J)=(4,2)$ gyromoment simulations lead to a smoothing of the distribution functions, reducing the sharp feature that appears around the thermal velocity (see figure 7, around $s_{\parallel a}=1$). This feature can be seen also in the lowest drive simulation in figure 7 where the $(P,J)=(10,5)$ gyromoment results present finer structures than the $(4,2)$ resolution. This smoothing effect confirms the hypothesis that the use of a reduced number of gyromoments yields an effective diffusion in the velocity space.

Figure 7. Time-averaged normalized ion velocity distribution function $|g_{v,i}(s_{\parallel,i},x_i)/ g_{v,i}(0,0)|$. The results from GENE (a,d) and from the gyromoment approach $(P,J)=(4,2)$ (b,e), $(P,J)=(10,5)$ (c,f), are presented for $\kappa _N=2.5$ (a–c) and $\kappa _N=1.6$ (d–f) keeping $\eta =0.25$ and $\nu =0$.

In conclusion, we remark that the gyromoment method shows its ability to simulate the Z-pinch nonlinear turbulent dynamics in the collisionless limit, which represents the most challenging regime for this approach. Valid results are obtained at large gradient drives with a velocity space represented by only 9 gyromoments per species, compared with the, approximately, $10^2$ velocity grid points used in GENE simulations. On the other hand, at the weakest gradient drive studied, convergence is obtained with a number of gyromoments approximately equal to the number of points used by GENE. In all cases, results obtained with a lower number of gyromoments still provide a reasonable prediction of the time-averaged level of transport.

4.1. Collisions and convergence

Adding collisions to our system helps significantly the convergence of the moment approach. Figure 8 shows the linear growth rates of the entropy mode for various Hermite–Laguerre bases, two collision frequency values, $\nu =0.01$ and $0.1$, and two gradient levels, $\kappa _N=1.6$ and $2.2$. This illustrates that convergence is obtained at high collisionality and high gradient levels with a low number of polynomials.

Figure 8. Convergence study of the entropy mode growth rate for $\kappa _N=2.2$ (a,b) and $\kappa _N=1.6$ (c,d) using the GK Sugama collision operator with $\nu = 0.01$ (a,c) and $\nu =0.1$ (b,d) for $\eta = 0.25$. The colour indicates the polynomial basis used: $(P,J)=(2,1)$ (cyan), $(P,J)=(4,2)$ (blue), $(P,J)=(6,3)$ (pink), $(P,J)=(8,4)$ (red) and $(P,J)=(10,5)$ (black).

Similarly, figure 9 shows the nonlinear transport level for the parameters of figure 8. One can observe, in particular, that the $(P,J)=(2,1)$ basis is sufficient at high collisionality and high gradient values. In the other cases, nonlinear simulations carried out with this reduced polynomial basis overestimate the level of transport when the linear growth rate is overestimated if evaluated with the same number of polynomials, and vice versa. Finally, the linear and nonlinear results presented in figures 8 and 9, respectively, demonstrate that the basis $(P,J)=(4,2)$ is sufficient to obtain accurate results in the parameter region of interest. Thus, the linear and nonlinear collisional simulations are performed using the polynomial basis $(P,J)=(4,2)$ in the following. As an indication of the computational cost of our nonlinear simulations, we notice that one Runge-Kutta 4 time-step for $(P,J)=(4,2)$ and $200\times 64$ spatial points is performed, on average, in $48$ ms (wall clock time) when run on one Marconi node, i.e. $2\times 24$-cores Intel Xeon 8160 (SkyLake) at 2.10 GHz.

Figure 9. Convergence study of the turbulent transport time traces for $\kappa _N=2.2$ (a,b) and $\kappa _N=1.6$ (c,d) using the GK Sugama collision operator with $\nu = 0.01$ (a,c) and $\nu =0.1$ (b,d). The colour indicates the polynomial basis used, $(P,J)=(2,1)$ (cyan), $(P,J)=(4,2)$ (blue): $(P,J)=(6,3)$ (pink), $(P,J)=(8,4)$ (red). The other parameters are $\eta = 0.25$ and $N_x=200$, $N_y=64$ for the spatial resolution.

4.2. Impact of collisions on the entropy mode and the Dimits shift

Figure 10 shows the impact of collisions on the entropy mode linear growth rate for the cases considered in § 3. Collisions stabilize the tail of the entropy mode present at high $k_y$ in the collisionless regime because of diffusion in phase space, as observed in Ricci et al. (Reference Ricci, Rogers, Dorland and Barnes2006b). This effect is recovered for both collision frequencies and by all the operators considered here, which also include GK effects that induce strong damping for $k_y\gtrsim 1$. At low $k_y$ and in the proximity of its peak value, one can observe that the growth rate is affected by collisions in different ways, depending on the collision model. On the one hand, large-scale fluctuations are destabilized by collisional effects in the case of the Dougherty and Sugama collision operators for $\kappa _N=2.0$ and $\kappa _N=2.5$. In this case, an increase of the growth rate at $k_y\sim 0.5$ is observed. This effect is similar to the one observed in instabilities that have a fluid nature, such as the drift waves, which are destabilized by resistivity (Goldston & Rutherford Reference Goldston and Rutherford1995). On the other hand, the collisional growth rate is smaller or close to the collisionless case when the Coulomb and Lorentz collision operators are used, for all values of $\kappa _N$ and $k_y$ considered.

Figure 10. Linear growth rate of the entropy mode for different collision models and comparison with the collisionless results (black) for two different collision frequencies, $\nu =0.1$ (a–c) and $\nu =0.01$ (d–f) and for three different drive values, $\kappa _N=1.6$ (a,d), $\kappa _N=2.0$ (b,e) and $\kappa _N=2.5$ (c,f), keeping $\eta =0.25$. The different lines denote the Dougherty (red), Sugama (blue), Lorentz (yellow) and Coulomb (green) operators used in the gyromoment approach with a $(4,2)$ Hermite–Laguerre basis.

We now turn to the nonlinear results that include finite collisionality, and we discuss two scans of simulations, for $\nu =0.1$ and $\nu =0.01$, where the drive value is varied from $\kappa _N=1.6$ to $\kappa _N=2.5$, being $\eta = 0.25$. The results are shown in figure 11, where they are compared with the collisionless limit, which shows a Dimits threshold value $\kappa _N\simeq 2$, similarly to Hallenbert & Plunk (Reference Hallenbert and Plunk2022), below which ZF suppress turbulence. We observe that the effect of collisions vanishes at large drive values, where the ZF do not play a crucial role, in agreement with the observations in Ricci et al. (Reference Ricci, Rogers and Dorland2006a). This suggests that the effect of collisions is mostly related to the ZF dynamics and, as we show later, through their damping and related weakening of the associated transport barrier.

Figure 11. Collisional saturated transport level for different collision operators at $\nu =0.01$ (a) and $\nu =0.1$ (b): Dougherty (red triangles), Sugama (blue squares), modified Sugama (light blue squares), Coulomb (green diamonds) and Lorentz (yellow triangles). The collisionless results are also reported (black stars) with the mixing length estimate $\varGamma _x^\infty \sim \gamma _{p}^2/k_p^3$ (dashed black line). In all cases, $\eta =0.25$.

When turbulence is fully developed, the amplitude of the fluctuations can be estimated considering a balance between the nonlinear saturating terms and the linear drive, $\partial _t \sim \boldsymbol v_{E\times B}\boldsymbol {\cdot } \boldsymbol {\nabla }$. This yields $\gamma \sim k^2 \phi$ and, thus, $\phi \sim \gamma /k^2$, considering the peak linear growth rate, $\gamma$, the associated wavenumber, $k$, and with the assumption of circular eddies ($k_x \sim k_y \sim k$). Using the Poisson equation, one observes that the particle density scales with the potential fluctuations, $n\sim \phi$, which leads to the estimate of the radial particle transport $\varGamma _x\sim \gamma ^2/k^3$. This scaling, based on the collisionless peak value of the entropy mode instability, is shown in figure 11, revealing that it captures well the dependence of $\varGamma _x$ at strong gradients, where the effect of ZF is weak. On the other hand, the reduction of the transport by the ZF cannot be captured by this mixing-length estimate at a low gradient value.

At medium and low levels of the driving gradient, where ZF are expected to play a role, according to the collisionless values, the different collision models lead to significantly different results, for both the $\nu = 0.01$ and $\nu =0.1$ cases. In particular, the Sugama and Dougherty tend to differ from the Lorentz and Coulomb operators. The difference cannot be explained solely in terms of linear growth rate since the Coulomb operator linear results differ from the Lorentz results at lower gradient values (see figure 10). In fact, the ZF quenching of the turbulence (Kobayashi & Rogers Reference Kobayashi and Rogers2012) has a strong dependence on the collision model.

Decreasing the collision frequency by a factor of ten, i.e. between $\nu =0.1$ and $\nu =0.01$ (see figure 11), reduces the gap between collisional and collisionless results in the majority of parameters and collision models studied. However, it is worth noting that, at a high gradient value, the transport does not approach the collisionless value monotonically with resistivity. This phenomenon is due to a combination of the tertiary instability affecting the ZF and the damping of turbulence due to collisions. In fact, Ricci et al. (Reference Ricci, Rogers and Dorland2006a) observed a non-monotonic dependence of transport to collisionality at large gradient values as well. We note that this feature does not depend on the chosen collision model.

The results obtained with the Dougherty operator appear to most closely approach the collisionless case, with Dougherty being the only operator that shows a Dimits shift at $\kappa _N\simeq 2.1$ for $\nu = 0.1$ and $\kappa _N \simeq 2.3$ for $\nu = 0.01$. This similarity can be explained by the simplicity of the Dougherty model, which is mainly composed of kinetic and spatial diffusion terms that are present, albeit at smaller amplitude and for numerical reasons, also in the collisionless case. Concerning the collisionless case, we expect that the slight reduction of transport at the lowest drive level is due to the reduced linear drive. When reducing the collisionality at $\kappa _N\leq 1.7$, we observe that transport is significantly reduced in comparison with the collisionless case. In this regime, the ZF are stable and yield small bursts of transport occurring over large time intervals, approximately $2000 L_B/c_s$. We report that increasing the polynomial basis to $(P,J)=(8,4)$ does not affect significantly the result obtained with the $(P,J)=(4,2)$ basis at these lower gradient values.

Similarly to the Dougherty operator, the Sugama operator yields a regime of suppressed transport at low gradient values and a regime of fully developed turbulence at large gradient values. However, at intermediate gradient levels, transport is remarkably larger in comparison with the collisionless results and Dougherty operator, with the oscillations between quiescent and turbulent periods (see figure 3) being replaced by fluctuations around a plateau value with persistent ZF structures. This feature, also observed with the Lorentz and Coulomb operators, can be explained by a ZF damping sufficiently strong to continuously allow fluctuations to grow in the ZF regions where the $\boldsymbol E\times \boldsymbol B$ velocity shear vanishes (Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). In the context of the predator–prey cycles, this case corresponds to an overlap of bursts. This effect is reduced when the collision frequency is decreased to $\nu =0.01$ where cyclic transport dynamics, previously identified by Kobayashi & Gürcan (Reference Kobayashi and Gürcan2015) and shown in figure 12, is obtained. The frequency of these bursts is directly related to the ZF damping rate due to the collision operator that dissipates the ZF structures (highlighted by the decreasing phase of the zonal energy, blue line of figure 12), and the primary instability growth rate (underlined by the slope of the increasing part of the non-zonal energy, red line in figure 12). The fact that the burst period is shorter with the Sugama operator than the Dougherty operator indicates that a strong ZF damping mechanism resides in the higher gyromoment coupling present in the Sugama operator.

Figure 12. Amplitude of the normalized spectral energy for zonal modes ($\sum _{k_y=0}|\phi _k|$; blue), non-zonal modes ($\sum _{k_y\neq 0}|\phi _k|$; red) and transport level ($\varGamma _x$; yellow) obtained for a simulation with the Sugama collision operator for $\kappa _N=1.6$, $\eta =0.25$ and $\nu = 0.01$, with $(P,J)=(4,2)$ gyromoments.

The reduction of the transport level with respect to the mixing length estimate at the lowest gradients is less pronounced with the Lorentz operator than with the Dougherty and Sugama operators. The difference between the Sugama and the Lorentz model is mainly due to the energy diffusion term contained in the field part of the former collision operator. In fact, the pitch-angle scattering Lorentz operator does not contain any energy diffusion term, while the Sugama model uses an ad hoc energy diffusion term in the field part of the collision operator (on the other hand, the spatial diffusion terms of the Lorentz and Coulomb operators coincide). Confirming the importance of having an accurate description of the energy diffusion, we note that tests at low drive values, where we modify the Sugama operator by zeroing out the ad hoc energy diffusion term (while also breaking the Sugama conservation properties), show a significant increase of the transport level (light blue squares in figure 11). It is worth noting that no bursts are observed in Lorentz simulations even in the low collisionality case.

The Coulomb collision operator simulations do not show remarkable differences in comparison with the Lorentz collision operator in the high collisionality case. Both operators maintain a high level of transport, even at low gradient values. It is worth noting that the Coulomb collision operator induces the largest level of transport than all other collision operators for almost every $\kappa _N$, which can be surprising since the related linear growth rate is smaller than the one yielded by the other collision operators (see figure 10). In particular, in the low collisionality case the Coulomb operator maintains a high transport level also with respect to the Lorentz operator.

Confirming that collisions regulate transport through the ZF damping, we now describe a detailed study of this mechanism, as induced by the different collision models. We consider the nonlinear collisionless saturated states for $\kappa _N = 1.6$, 2.0, and 2.5 (see figure 5) as initial conditions for a set of simulations that use different collision models. We isolate the damping effect by removing the entropy mode drive, $\kappa _N=0$, and we use a $(P,J)=(4,2)$ gyromoment set with $\nu =0.1$. We let the system evolve and follow the damping of the ZF profile. The results of this numerical experiment can be first observed qualitatively in figure 13 where the averaged radial ZF profile, $\langle \partial _x\phi \rangle _y$, is plotted as a function of time for each collision operator considered. Figure 13 reveals that the effect of collisions on the ZF profile is highly dependent on the operator model. The Dougherty model does not significantly affect the ZF structure, while the Sugama operator leads to their damping. The Lorentz operator filters the initial ZF structure, decreasing the amplitude of short-wavelength ZF, while a long-wavelength mode survives. Finally, the Coulomb operator strongly damps the ZF at all wavelengths. Thus, confirming our hypothesis that the different ZF damping is responsible for the different level of transport, the smallest transport values observed in figure 11 correspond to the operators that allow the smallest scale of the ZF structure to survive.

Figure 13. Time evolution of the $y$-averaged ZF profile, $\langle \partial _x\phi \rangle _y$, for the Dougherty (a), Sugama (b), Lorentz (c) and Coulomb (d) collision operators, using the saturated state of the collisionless simulation at $t_0=5000$ for $\kappa _N=1.6$, $\nu =0.1$ and $\eta =0.25$ as initial conditions.

As a further confirmation and a more quantitative analysis of the results shown in figure 13, we define the normalized ZF energy, i.e.

(4.1)\begin{equation} A^2_{ZF}(t) =\frac{\displaystyle \int {{\rm d} x} \langle\phi\rangle_{y}^2(t)}{\displaystyle \int {{\rm d} x} \langle\phi\rangle_{y}^2(0)}, \end{equation}

and study its time evolution for each collision operator in figure 14. As initial conditions, we consider the ZF obtained in the $\kappa _N=1.6$, $\kappa _N=2.0$ and $\kappa _N=2.5$ collisionless simulations. Focusing on the damping at early times, $\partial _t A^2_{ZF}|_{t=0}$ (the growth of the linear instability alters the ZF damping at time scales $1/\gamma \sim 10$), this analysis unveils a clear difference between Dougherty and Sugama operators, while these operators provide very similar linear growth rates. We also observe that the Lorentz and Coulomb operators yield similar damping, corresponding to a similar transport level in the nonlinear simulations. Thus, we can deduce that, unlike the linear growth rate, the saturated transport level is directly related to the ZF damping rate.

Figure 14. Time evolution of $A_{ZF}^2$ (see (4.1)) for the Dougherty (red), Sugama (blue), Lorentz (yellow) and Coulomb (green) collision operators used in the gyromoments approach with a $(4,2)$ Hermite–Laguerre basis. The ZF initial conditions are the ones obtained from the $\kappa _N=1.6$ (a), $\kappa _N=2.0$ (b) and $\kappa _N=2.5$ (c) collisionless simulations. In all cases $\eta =0.25$.