Gyrokinetic simulations of plasma turbulence in a Z-pinch using a moment-based approach and advanced collision operators

Abstract

The first nonlinear gyrokinetic simulations obtained using a moment approach based on the Hermite–Laguerre decomposition of the distribution function are presented, implementing advanced models for the collision operator. Turbulence in a two-dimensional Z-pinch is considered within a flux-tube configuration. In the collisionless regime, our gyromoment approach shows very good agreement with nonlinear simulations carried out with the continuum gyrokinetic code GENE, even with fewer gyromoments than required for the convergence of the linear growth rate. By using advanced linear collision operators, the role of collisions in setting the level of turbulent transport is then analysed. The choice of collision operator model is shown to have a crucial impact when turbulence is quenched by the presence of zonal flows. The convergence properties of the gyromoment approach improve when collisions are included.

1. Introduction

The understanding of the dynamics in the tokamak boundary, the region that encompasses the edge and scrape-off-layer (SOL), is crucial to predict the performance of future tokamak devices. Although gyrokinetic (GK) models provide accurate simulations of the plasma dynamics in core conditions, they come at a high computational cost. Alternatively, fluid approaches are computationally less expensive but, being limited to evolving a finite number of moments of the distribution function, their results are questionable when the hypothesis behind the closure (such as the high collisionality closure) do not hold. To overcome the limitations of current models, Frei, Jorge & Ricci (2020) propose an extension of the SOL drift-kinetic model presented in Jorge, Ricci & Loureiro (2017) and develop a GK model based on the projection of the velocity space dependence of the distribution function onto a Hermite–Laguerre polynomial basis. Extending previous full-F gyrofluid models (Strintzi, Scott & Brizard 2005; Madsen 2013; Held, Wiesenberger & Kendl 2020), for instance in the number of evolved gyromoments, this set of fluid equations converges to the description of the evolution of the distribution function provided by the full-F GK Boltzmann equation, as the number of moments increases. In the Hermite–Laguerre framework, advanced collision operators such as the full nonlinear Coulomb collision operator (Jorge, Frei & Ricci 2019) as well as linearized ones (Frei et al. 2021), can be used to model collisional effects. In the $\delta f$ limit of the full-F model presented by Frei et al. (2020), the Hermite–Laguerre decomposition can be interpreted as an extension of the previous gyrofluid model (Brizard 1992; Hammett, Dorland & Perkins 1992; Beer, Cowley & Hammett 1995; Snyder & Hammett 2001; Scott 2005) to an arbitrary number of moments. This approach, also pursued by Mandell, Dorland & Landreman (2018a) with the GX code (Mandell, Dorland & Landreman 2018b; Mandell et al. 2022), yields an infinite set of fluid equations for the basis coefficients, the gyromoments, which describe the deviations of the distribution function from a Maxwellian distribution. The efficiency of the $\delta f$ Hermite–Laguerre gyromoment approach is demonstrated by Frei, Hoffmann & Ricci (2022a) focusing on the linear properties of the ion temperature gradient (ITG) instability in the slab limit, as well as in a flux-tube geometry (Frei et al. 2022b), including the use of the linearized GK Landau form of the Fokker–Planck collision operator. These works demonstrate the improvement of convergence properties of the gyromoment method with collisions, i.e. when deviations from a Maxwellian distribution function are reduced. In addition, even in the collisionless case, it is shown that the number of gyromoments needed for linear convergence is less than the number of grid points necessary for convergence in the state-of-the-art continuum GK code GENE (Jenko, Dorland & Kotschenreuther 2000).

Here, nonlinear simulations are presented for the first time using a Hermite–Laguerre gyromoment approach. We consider a local Z-pinch geometry which is characterized by a cylindrically symmetric plasma confined by a purely azimuthal, radially dependent, magnetic field with equilibrium radial gradients in temperature and density. In the presence of a background density gradient, an entropy mode (Ricci et al. 2006b) develops in the Z-pinch that can be modelled by using a local $\delta f$ GK approach with a kinetic treatment of the electrons. This mode develops perpendicularly to the magnetic field and persists in the $k_\parallel =0$ limit, where $k_\parallel$ denotes the perturbation wavenumber along the magnetic field, allowing simulations to be performed for a limited computational cost. While the Z-pinch geometry is considerably simpler than, e.g. the one in a tokamak (for instance, it does not have magnetic shear nor toroidal effects such as particle trapping), it still allows for the study of complex nonlinear phenomena, such as the emergence of zonal flows (ZF) (Fujisawa et al. 2004; Diamond et al. 2005) that lead to the Dimits shift (Dimits et al. 2000), which role continues to challenge our understanding of tokamak physics.

The first nonlinear simulations in a Z-pinch, presented by Ricci, Rogers & Dorland (2006a), study the level of transport induced by the entropy mode as a function of the density gradient, showing that ZF can regulate the level of turbulent transport. However, the effect of ZF can be reduced either as the result of collisions, modelled in Ricci et al. (2006a) through a drift-kinetic (DK) Lorentz operator, or by a tertiary Kelvin–Helmholtz instability (KHI), destabilized in scenarios characterized by a sufficiently large density gradient drive. These results are confirmed by Kobayashi & Rogers (2012) using a GK single-species collision operator described in Abel et al. (2008) and Barnes et al. (2009).

At low-density gradient drive, i.e. under the tertiary KHI instability threshold, transport regimes characterized by bursts rising from the competition between ZF collisional damping and quenching of the primary instability are identified and modelled with a predator–prey cycle by Kobayashi, Gürcan & Diamond (2015). In order to explore the mechanisms behind the ZF formation and damping, Ivanov et al. (2020) use a fluid-diffusive collision operator obtained by integration of the linearized Coulomb collision operator and derive a three-field, two-dimensional fluid model directly from the GK equation in a Z-pinch geometry, later extended to three dimensions (Ivanov, Schekochihin & Dorland 2022). This model includes first-order finite Larmor radius (FLR) effects in the long-wavelength, cold-ion limit and allows exploration of the ZF dynamics within an analytical framework. The simulations show good qualitative agreement with modified Hasegawa–Wakatani simulations (Qi, Majda & Cerfon 2020). Similarly, Hallenbert & Plunk (2021) derive a fluid model in a Z-pinch geometry in the collisionless limit, including second-order FLR effects. This allows the numerical prediction of the Dimits threshold, i.e. the gradient level below which transport is strongly reduced by the presence of ZF (Hallenbert & Plunk 2022). The prediction is confirmed by comparison with GENE simulations.

The present paper reports on the first nonlinear GK simulations carried out with the Hermite–Laguerre gyromoment approach using advanced collision operators (Hoffmann & Frei 2020). These simulations include nonlinear $\boldsymbol E \times \boldsymbol B$ advection, FLR effects of arbitrary order, kinetic electrons and, leveraging the work in Frei et al. (2021), a set of advanced linear GK collision operators. These operators include the single-species Dougherty model (Dougherty 1964), the multi-species Sugama model (Sugama, Watanabe & Nunami 2009), the single-species pitch-angle scattering operator with a restoring momentum term, denoted as the Lorentz operator (Helander & Sigmar 2002), and the Landau form of the multi-species Fokker–Planck model that we denote as the Coulomb operator (Rosenbluth & Longmire 1957; Hazeltine & Meiss 2003). We consider a local $\delta f$ flux-tube approach that separates equilibrium and fluctuating quantities, assuming constant equilibrium gradients across the domain. By imposing $k_\parallel =0$, we evolve the turbulent dynamics on a perpendicular plane. This set-up provides an ideal framework to compare the gyromoment model with a continuum code in a nonlinear turbulent regime, and to study the effect of advanced linearized collision models in ZF-dominated systems.

Our results demonstrate, first, the ability of the gyromoment approach to retrieve linear and nonlinear collisionless results obtained with the GK continuum code GENE. In particular, we observe that the number of gyromoments needed for convergence increases while approaching the linear marginal stability conditions, and that underresolved collisionless simulations present predator–prey cycles, typically observed in collisional GK simulations (Kobayashi et al. 2015) and fluid-reduced models (Qi et al. 2020). The same dynamics is observed when increasing significantly the numerical dissipation acting on the velocity space in GENE. Secondly, we present a set of simulations at different instability drives in the collisionless limit and in the presence of collisions, which are modelled using the Dougherty, Sugama, Lorentz and Coulomb collision operators. The particle flux reveals a Dimits threshold in the collisionless limit. For gradient levels above the Dimits threshold and at finite collisionality, we observe negligible differences between the different collision operators. Shear flow stabilization effects are negligible and turbulence is fully developed. The transport is well approximated by a mixing-length argument, $\varGamma _x\sim \gamma ^2/k^3$ (Ricci et al. 2006a), where $\varGamma _x$ is the saturated particle transport level along the radial direction, while $\gamma$ and $k$ are the peak linear growth rate and wavelength of the entropy mode, respectively. Below the Dimits threshold, turbulence is quenched by ZF, which may be damped by collisions, and the choice of collision model affects significantly the transport level. A study of the ZF collisional damping provides and explanation for the differences observed between the collision operators.

The paper is organized as follows. In § 2, we briefly describe the nonlinear GK model in Z-pinch geometry and develop the gyromoment approach in this configuration. Section 3 presents linear and nonlinear benchmarks of the gyromoment approach with GENE in the collisionless limit. The dependence of the transport level with the instability drive and the role of collisions is investigated in § 4. The conclusions follow in § 5. In Appendix A, we show that a gyrofluid model as well as an extended Hasegawa–Wakatani model can be obtained by properly truncating the gyromoment equation hierarchy.

2. Gyrokinetic model of a Z-pinch configuration based on the gyromoment model

In this section we present, first, the GK model in the Z-pinch geometry, considering the local $\delta f$ flux-tube limit. Second, we project the Z-pinch GK equation on a Hermite–Laguerre polynomial basis in velocity space, thus obtaining an infinite set of two-dimensional equations for the gyromoments, which we denote as the gyromoment equation hierarchy. Finally, we present the numerical implementation of this hierarchy of equations.

2.1. Gyrokinetic model in a Z-pinch configuration

We consider the GK approach (Catto 1978; Frieman & Chen 1982; Hazeltine & Meiss 2003) 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 2003).

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 2007)

(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. 2006a; Ivanov et al. 2020), 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 2017; Frei et al. 2020), 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 2014). 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 1967; Grant & Feix 1967; Joyce, Knorr & Meier 1971; Gibelli & Shizgal 2006; Parker & Dellar 2015). 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. 2017), their main application is in expressing collision models in a spectral framework (Brunner, Valeo & Krommes 2000; Belli & Candy 2012).

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. 2020). 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. 2018a, 2022).

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

(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 (1992) can be retrieved by properly truncating the collisionless gyromoment hierarchy. We also show how an extended Hasegawa–Mima model (Hasegawa & Mima 1978; Dewhurst, Hnat & Dendy 2009) 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. 2018a).

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 2002; Hazeltine & Meiss 2003), 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. 2021)

(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. 2009) is a multi-species generalization of the Abel operator (Abel et al. 2008). 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. (2006a) 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 1958) 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. (2021). 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 (2011) and Loureiro et al. (2016), 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 1991). 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 2005) 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 1971).

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. (2021). 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. (2021) for more details). We ensured the convergence of our matrix evaluation by analysing the eigenvalue spectrum and the matrix symmetry.

3. Collisionless limit and comparison with the GENE code

In the present section, we analyse the results of the gyromoment simulations in the collisionless limit and demonstrate the ability of this approach to retrieve the results of the continuum GK code GENE in a collisionless two-dimensional Z-pinch configuration considering an equilibrium plasma with $\tau =1$, and a realistic electron–proton mass ratio, $\sigma = \sqrt {m_e/m_i}= 0.023$. GENE simulations are set up by closely following Hallenbert & Plunk (2022).

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. (2006b). 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. (2006b), 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. 2006b). 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. 2022a) and the toroidal ITG instability (Frei et al. 2022b). 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 2022). 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 (2022) are retrieved. It is worth noting that Hallenbert & Plunk (2022) 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 2010). 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 (2022). 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 2005), 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. (2006a) and Hallenbert & Plunk (2022). 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 1958). 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. Collisional turbulent transport

Building on the benchmark of our gyromoment solver with the GENE code in the collisionless limit, we now study the gyromoment method at finite collisionality, in particular $\nu =0.1$ and $\nu = 0.01$. These values encompass the typical collision rate in the core of a tokamak device (e.g. the collision frequency estimate in the DIII-D cyclone base case corresponds to $\nu \sim 0.05$ in our normalized units Lin et al. 1999). We first present the impact of collisions on the convergence of the Hermite–Laguerre basis using the Sugama collision operator. Then, we investigate the properties of turbulence in a Z-pinch, as obtained by using different linear collision operators.

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. (2006b). 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 1995). 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 (2022), 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. (2006a). 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 2012) 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. (2006a) 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. 2020). 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 (2015) 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$.

5. Conclusions

In the present paper, the first nonlinear GK simulations carried out using a gyromoment approach and including advanced collision models are presented. By implementing the moment hierarchy in (2.14), turbulence in a two-dimensional Z-pinch geometry is studied.

We first present a benchmark with the continuum GK code GENE that demonstrates the ability of the gyromoment approach to simulate accurately the nonlinear evolution of the entropy mode, even in the collisionless limit. We show that the convergence behaviour of the nonlinear results follows the same trend as the linear ones, i.e. convergence properties improve with the increase of the gradient strength. However, accurate nonlinear results require only that the linear growth rate of the modes developing at large scales are accurately resolved.

We then extend the nonlinear results, adding collisions with the use of four different collision operator models. We observe that the gyromoment simulations converge with a lower number of gyromoments than in the collisionless case. With a Dimits threshold identified around $\kappa _N\sim 2$ in the collisionless case, the influence of collisions on the transport level becomes particularly evident for $\kappa _N<2$. This confirms previous studies (Lin et al. 1999; Ricci et al. 2006a; Ricci, Rogers & Dorland 2010) pointing that collisional effects are mainly related to the dynamics of the ZF. Our results highlight the disagreement between Dougherty, Sugama, Lorentz and Coulomb GK collision models, in the linear growth rate and, even more, in the level of nonlinear transport. We show that the analysis of the linear results is not sufficient to predict the difference observed in the saturated transport level. However, we observe a direct link between ZF damping and transport level, which could be used to develop a reduced model of the transport level in a future work. By demonstrating for the first time that the transport level in ZF dominated regime is highly dependent on the collision model in use, we point out that the choice of the collision operator should be properly considered in GK turbulence simulations where ZF are present.

In a more general context, the present study is a first step towards the nonlinear simulation of the tokamak boundary based on the use of the gyromoment approach. Our plan is to consider nonlinear simulations in the $s-\alpha$ flux-tube geometry as a next step, expanding the linear study presented in Frei et al. (2022b). Using the cyclone base case as a reference (Dimits et al. 2000), we will have a benchmark for evaluating the performance of the gyromoment approach. In particular, the ability of the Hermite–Laguerre approach to accurately resolve nonlinear trapped-particle dynamics remains an open area of research. We will then turn to the simulation of the tokamak boundary. According to Frei et al. (2022b), we anticipate that the high pressure gradient and level of collisionality present in the tokamak edge will improve the convergence of the nonlinear gyromoment hierarchy with respect to core conditions. However, the $\delta f$ assumption will need to be relaxed when simulating the SOL. The convergence behaviour of a Hermite–Laguerre moment approach in a full-F nonlinear framework and the influence of its closure model remains a central question.