2.1. Gyrokinetic Vlasov–Maxwell system

The GENE-3D code is a Eulerian code that solves the gyrokinetic Vlasov equation, coupled to Maxwell's equations (e.g. Brizard & Hahm Reference Brizard and Hahm2007) on a five-dimensional grid in phase space, using two velocity dimensions – the velocity parallel to the equilibrium magnetic field $v_{||}$ and the magnetic moment $\mu =mv_{\perp }^2/2B_0$ – and three spatial dimensions $(x,y,z)$, representing a field-aligned coordinate system (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995). Furthermore, it employs the so-called ‘$\delta f$’ splitting (e.g. Garbet et al. (Reference Garbet, Idomura, Villard and Watanabe2010), and references therein). In this approach, one assumes that the full gyrocentre distribution function $F_{{\sigma }}$ of species $\sigma$ can be split into a stationary background distribution function $F_{0,\sigma }$, describing the plasma at thermal equilibrium, and a time-dependent, first-order perturbation $F_{1,\sigma }$:

(2.1)\begin{equation} F_{\sigma}=F_{0,\sigma}+F_{1,\sigma},\quad \mathrm{with} \ ||F_{1,\sigma}||/||F_{0,\sigma}|| \sim O \left( \epsilon_{\delta} \right)\ll 1, \end{equation}

where $\epsilon _{\delta }$ is a smallness parameter describing the scale relation between background and perturbed quantities. The GENE-3D code assumes that the background distribution is given by a local Maxwellian:

(2.2)\begin{equation} F_{0,\sigma}=F_{M,\sigma}({\boldsymbol{x}},v_{||},\mu)=\frac{n_{0,\sigma}({\boldsymbol{x}})}{({\rm \pi})^{{3}/{2}} v_{\textrm{th},\sigma}({\boldsymbol{x}})^{3}} \exp \left( -\frac{m_{\sigma}v_{||}^2/2 +\mu B_{0}({\boldsymbol{x}})}{T_{0,\sigma}({\boldsymbol{x}})}\right). \end{equation}

Here, $m_{\sigma }$, $n_{0,\sigma }$ and $T_{0,\sigma }$ are the mass, equilibrium background density and temperature of particle species $\sigma$ at position ${\boldsymbol {x}}$, respectively. Furthermore, $B_0$ gives the strength of the equilibrium magnetic field, $v_{||}$ and $\mu$ denote the velocity parallel to the magnetic field and the magnetic moment, respectively, and $v_{\textrm {th},\sigma }=\sqrt {2T_{0,\sigma }({\boldsymbol {x}})/m_{\sigma }}$ is the thermal velocity of particle species $\sigma$. Under this assumption, the gyrokinetic equation for particle species $\sigma$, including a perturbation of the parallel vector potential, can be written as

(2.3)\begin{align} \frac{\partial F_{1,\sigma}}{\partial t}& ={-} \left[ v_{||} \boldsymbol{b}_0 +\left( \boldsymbol{v}_{\chi}+ \boldsymbol{v}_{\boldsymbol{\nabla} B}+\boldsymbol{v}_{c} \right) \right] \boldsymbol{\cdot} \boldsymbol{\nabla} F_{1,\sigma} +\frac{\mu}{m_{\sigma}} \boldsymbol{b}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} B_0 \frac{\partial F_{1,\sigma}}{\partial v_{||}} \nonumber\\ & \quad- \boldsymbol{v}_{\chi} \boldsymbol{\cdot} \left[ \boldsymbol{\nabla} \ln (n_{0,\sigma})+ \boldsymbol{\nabla} \ln (T_{0,\sigma}) \left( \frac{m_{\sigma} v_{||}^2/2 +\mu B_0}{T_{0,\sigma}}- \frac{3}{2} \right)\right] F_{M,\sigma}\nonumber\\ & \quad-\frac{q_{\sigma} F_{M,\sigma}}{T_{0,\sigma}} \left[ v_{||} \boldsymbol{b}_0 + \left( \boldsymbol{v}_{\chi}+\boldsymbol{v}_{\boldsymbol{\nabla} B}+\boldsymbol{v}_{c} \right) \right] \boldsymbol{\cdot} \boldsymbol{\nabla} \mathcal{G}\{ \phi_1 \} - \frac{q_{\sigma} v_{||}}{c} \frac{F_{M,\sigma}}{T_{0,\sigma}} \frac{\partial \mathcal{G}\{A_{1,||}\}}{\partial t}\nonumber\\ & \quad - \left(\boldsymbol{v}_{\boldsymbol{\nabla} B}+\boldsymbol{v}_{c} \right) \boldsymbol{\cdot} \left[ \boldsymbol{\nabla} \ln (n_{0,\sigma})+ \boldsymbol{\nabla} \ln (T_{0,\sigma}) \left( \frac{m_{\sigma} v_{||}^2/2 +\mu B_0}{T_{0,\sigma}}-\frac{3}{2} \right)\right] F_{M,\sigma}\nonumber\\ & \quad - \boldsymbol{v}_{E_0}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{1,\sigma}+C[F_{1,\sigma}], \end{align}

where $\chi = \phi _1 - v_{||} A_{1,||}/c$ is called the gyrokinetic potential. The operator $\mathcal {G}\{ . \}$ defines the forward-gyroaverage operator

(2.4)\begin{equation} \mathcal{G}\{u\}(\boldsymbol{X})=\frac{1}{2{\rm \pi}}\, \int\limits_0^{2{\rm \pi}}\, u(\boldsymbol{X}+ \boldsymbol{r}(\vartheta))\, \textrm{d}\vartheta, \end{equation}

where $\boldsymbol {X}$ denotes the gyrocentre position and $\boldsymbol {r}$ is the gyroradius vector orthogonal to the local magnetic field. The drift velocities appearing in (2.3) are defined as follows:

(2.5)\begin{equation} \left.\begin{gathered} \boldsymbol{v}_{E_0}\equiv\frac{c}{B_0^2} \boldsymbol{B}_0\times \boldsymbol{\nabla} \phi_0, \quad \boldsymbol{v}_{\chi}\equiv\frac{c}{B_0^2} \boldsymbol{B}_0 \times \boldsymbol{\nabla} \mathcal{G}\{ \chi\},\\ \boldsymbol{v}_{\boldsymbol{\nabla} B}\equiv \frac{\mu c}{q_{\sigma} B_0^2} \boldsymbol{B}_0 \times \boldsymbol{\nabla} B_0, \quad \boldsymbol{v}_{c}\equiv\frac{v_{||}^2}{\varOmega_{\sigma}}\left( \boldsymbol{b}_0 \times \left( \boldsymbol{\nabla} \ln(B_0) + \frac{\beta}{2} \boldsymbol{\nabla} \ln(p_0) \right) \right), \end{gathered}\right\} \end{equation}

where $\varOmega _{\sigma } =(q_{\sigma } B_0)/(m_{\sigma } c)$ is the gyrofrequency of species $\sigma$, $\boldsymbol {b}_0=\boldsymbol {B}_0/B_0$ is the unit vector in the direction of the equilibrium magnetic field $\boldsymbol {B}_0$, $p_0$ is the equilibrium thermodynamic pressure, $\beta =8 {\rm \pi}p_0/B_0^2$ is the ratio between plasma pressure and magnetic field strength and $c$ is the speed of light.

At the current stage, GENE-3D uses a linearised Landau–Boltzmann collision operator for $C[F_{1,\sigma }]$. Furthermore, the interplay between neoclassical and gyrokinetic physics is represented by the term coupling the curvature and $\boldsymbol {\nabla } B$ drifts to the background Maxwellian (Oberparleiter et al. Reference Oberparleiter, Jenko, Told, Doerk and Görler2016), and the effects of a long-wavelength, radial electric field on the system (Helander & Simakov Reference Helander and Simakov2008; Mishchenko & Kleiber Reference Mishchenko and Kleiber2012) are represented by the term proportional to $\boldsymbol {v}_{E_0}$ in (2.3). Though they are incorporated in the main algorithm, these effects are neglected in the simulations presented in this paper for simplicity.

The system is closed by the gyrokinetic field equations. Poisson's equation, providing the electrostatic potential, reads

(2.6)\begin{equation} \nabla^2 \phi_1(\boldsymbol{x},t) \approx \boldsymbol{\nabla}_{{\perp}}^2 \phi_1(\boldsymbol{x},t) ={-}4{\rm \pi}\, \sum_{\sigma}\, q_{\sigma} n_{1,\sigma}(\boldsymbol{x},t), \end{equation}

where only the spatial derivatives perpendicular to the magnetic field are retained, which is in accordance with the gyrokinetic ordering. In order to calculate the density perturbation at the particle position $\boldsymbol {x}$ from the gyrocentre position $\boldsymbol {X}$, GENE-3D employs a first-order pull-back transformation (e.g. Brizard & Hahm Reference Brizard and Hahm2007):

(2.7)\begin{equation} n_{1,\sigma}(\boldsymbol{x},t)=\int \mathcal{K}\{ F_{1,\sigma} \} -\left( \frac{q_{\sigma} F_{M,\sigma}}{T_{0,\sigma}} \phi_1 -\mathcal{K}\left\{ \frac{q_{\sigma} F_{M,\sigma}}{T_{ 0,\sigma}} \mathcal{G}\left\{\phi_1\right\} \right\} \right) \textrm{d}^3 v. \end{equation}

The operator $\mathcal {K}\{ . \}$ denotes the backward-gyroaverage operator, which is defined here as

(2.8)\begin{equation} \mathcal{K}\{u\}(\boldsymbol{x}) = \frac{1}{2{\rm \pi}}\, \int\limits_0^{2{\rm \pi}} u(\boldsymbol{x}- \boldsymbol{r}(\vartheta))\, \textrm{d}\vartheta. \end{equation}

Assuming a quasi-neutral limit ($\boldsymbol {\nabla }_{\perp }^2 \phi _1 \approx 0$), Poisson's equation can be rewritten as

(2.9)\begin{equation} \sum_{\sigma}\, q_{\sigma}^2 \, \int \left( \frac{F_{M,\sigma}}{T_{0,\sigma}} \phi_1 -\mathcal{K}\left\{ \frac{F_{M,\sigma}}{T_{0,\sigma}} \mathcal{G}\left\{\phi_1\right\} \right\} \right) \textrm{d}^3 v =\sum_{\sigma}\, q_{\sigma} \int\, \mathcal{K}\{ F_{1,\sigma} \}\, \textrm{d}^3v. \end{equation}

At the current stage, magnetic compressibility, associated with a parallel perturbation of the magnetic field, is neglected, so that only the parallel vector potential $A_{1,||}$ is retained in GENE-3D, which is determined by the parallel component of Ampere's law:

(2.10)\begin{equation} \boldsymbol{\nabla}_{\perp}^2\, A_{1,||} ={-}\frac{4 {\rm \pi}}{c} \sum_{\sigma}\, j_{1,||,\sigma} ={-}\frac{4 {\rm \pi}}{c} \sum_{\sigma}\,\int v_{||} \mathcal{K}\{ F_{1,\sigma} \}\, \textrm{d}^3v, \end{equation}

self-consistently with the gyrokinetic equation and Poisson's equation.

2.2. Electromagnetic model

Although (2.3), (2.9) and (2.10) are an analytically closed system, particular care has to be taken of the term containing the partial time derivative of $A_{1,||}$ in (2.3). The approach used in GENE for a long time was to introduce a new distribution function $g_{\sigma }=F_{1,\sigma }+(q_{\sigma } v_{||} F_{M,\sigma } A_{1,||})/( T_{0,\sigma } c )$ and solving the corresponding gyrokinetic system for $g_{\sigma }$ instead of $F_{1,\sigma }$ (Görler et al. Reference Görler, Lapillonne, Brunner, Dannert, Jenko, Merz and Told2011). However, it was shown in Crandall (Reference Crandall2019), that such a scheme becomes numerically unstable in nonlinear simulations at high plasma-$\beta$, although being stable linearly. In the same publication, an alternative scheme was proposed, similar to the one presented in Reynders (Reference Reynders1993), which has been included in the standard version of GENE recently. This scheme is also used in GENE-3D and is presented in the following. By defining the parallel inductive electric field as

(2.11)\begin{equation} E_{||}^\textrm{ind}\equiv{-}\frac{1}{c}\frac{\partial A_{1,||}}{\partial t}, \end{equation}

one can rewrite the gyrokinetic equation in the form

(2.12)\begin{equation} \frac{\partial F_{1,\sigma}}{\partial t} = R_{\sigma} +q_{\sigma} v_{||} \frac{F_{M,\sigma}}{T_{0,\sigma}}\mathcal{G}\left\{ E_{||}^{\textrm{ind}} \right\}, \end{equation}

where the term $R_{\sigma }$ contains all the contributions to the right-hand side of the gyrokinetic equation (2.3), except the one accounting for the parallel inductive electric field, and reads

(2.13)\begin{align} R_{\sigma}& ={-}\left[ v_{||} \boldsymbol{b}_0 +\left(\boldsymbol{v}_{\boldsymbol{\nabla} B}+ \boldsymbol{v}_{c} \right) \right] \boldsymbol{\cdot} \left( \boldsymbol{\nabla} F_{1,\sigma} + \frac{q_{\sigma} F_{M,\sigma}}{T_{0,\sigma}}\boldsymbol{\nabla} \mathcal{G}\{ \phi_1 \}\right)\nonumber\\ & \quad - \boldsymbol{v}_{\chi} \boldsymbol{\cdot} \left(\boldsymbol{\nabla} F_{1,\sigma} + \frac{q_{\sigma} F_{M,\sigma}}{T_{0,\sigma}}\boldsymbol{\nabla} \mathcal{G}\{ \phi_1 \} \right) +\frac{\mu}{m_{\sigma}} \boldsymbol{b}_0\boldsymbol{\cdot} \boldsymbol{\nabla} B_0 \frac{\partial F_{1,\sigma}}{\partial v_{||}}\nonumber\\ & \quad -\boldsymbol{v}_{\chi} \cdot \left[ \boldsymbol{\nabla} \ln (n_{0,\sigma})+ \boldsymbol{\nabla} \ln (T_{0,\sigma}) \left( \frac{m_{\sigma} v_{||}^2/2 +\mu B_0}{T_{0,\sigma}}-\frac{3}{2} \right)\right] F_{M,\sigma}\nonumber\\ & \quad - \left(\boldsymbol{v}_{\boldsymbol{\nabla} B}+\boldsymbol{v}_{c} \right) \boldsymbol{\cdot} \left[ \boldsymbol{\nabla} \ln (n_{0,\sigma})+ \boldsymbol{\nabla} \ln (T_{ 0,\sigma}) \left( \frac{m_{\sigma} v_{||}^2/2 +\mu B_0}{T_{0,\sigma}}-\frac{3}{2} \right)\right] F_{M,\sigma}\nonumber\\ & \quad - \boldsymbol{v}_{E_0}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{1,\sigma}+C[F_{1,\sigma}]. \end{align}

The next step is to take the partial time derivative of Ampere's law (2.10) and multiply by $(-1/c)$:

(2.14)\begin{equation} -\frac{1}{c}\boldsymbol{\nabla}_{{\perp}}^2\, \frac{\partial A_{1,||}}{\partial t} = \boldsymbol{\nabla}_{{\perp}}^2 E_{||}^{\textrm{ind}} = \frac{4 {\rm \pi}}{c^2} \sum_{\sigma}\,q_{\sigma}\,\int v_{ ||} \mathcal{K}\left\{ \frac{\partial F_{1,\sigma}}{\partial t} \right\} \textrm{d}^3v. \end{equation}

Inserting definition (2.12) into (2.14) eliminates the explicit coupling of the equation to the particle distribution function $F_{1,\sigma }$. Therefore, one obtains an equation that treats the inductive electric field as an independent electromagnetic field:

(2.15)\begin{align} \boldsymbol{\nabla}_{{\perp}}^2 E_{||}^{\textrm{ind}} - \frac{4 {\rm \pi}}{c^2} \sum_{\sigma}\,q_{\sigma}^2\,\int v_{||}^2\, \mathcal{K} \left\{ \frac{F_{M,\sigma}}{T_{0,\sigma}}\mathcal{G}\left\{ E_{||}^{\textrm{ind}} \right\} \right\} \textrm{d}^3v = \frac{4 {\rm \pi}}{c^2} \sum_{\sigma}\,q_{\sigma}\,\int v_{||} \mathcal{K}\{R_{\sigma} \}\, \textrm{d}^3v. \end{align}

Using (2.9), (2.10), (2.12), (2.13) and (2.15), the evolution of the distribution function $F_{1,\sigma }$ can be calculated by treating time as a discrete variable and using a numerical scheme to solve (2.12) in the above system approximately. For this, GENE-3D provides several explicit Runge–Kutta schemes, with a fourth-order-accurate Runge–Kutta scheme (RK4) being the default option.

3.1. Linear simulations with varying plasma-$\beta$

This subsection presents a linear verification study between GENE-3D and GENE over a large range of plasma pressures. For this, linear growth rates and mode frequencies of the most unstable mode are calculated varying the plasma-$\beta$ and are compared with the results presented in Görler et al. (Reference Görler, Tronko, Hornsby, Bottino, Kleiber, Norscini, Grandgirard, Jenko and Sonnendrücker2016).

The geometry is a tokamak with circular, concentric flux surfaces, with a magnetic field strength of $B_{0}=2.0\ \mathrm {T}$ on axis, which is called $B_{\textrm {ref}}$ in the following, a major radius of $R_0=1.67\ \mathrm {m}$, an inverse aspect ratio of $a/R_0=0.36$ and a safety factor profile

(3.1)\begin{equation} q(x)=0.86-0.16(x/a)+2.52(x/a)^2. \end{equation}

Furthermore, density and temperature profiles are defined according to

(3.2)\begin{equation} \frac{A}{A(x_0)}=\exp{\left[ -\kappa_{A} w_{A}\frac{a}{R_0} \tanh{\left( \frac{x-x_0}{w_{A} a}\right)} \right]} \end{equation}

for both electrons and ions, where $A=(n,T_{i},T_{e})$. The corresponding normalised gradients are

(3.3)\begin{equation} \frac{R_0}{L_{A}}={-}R_0 \frac{\partial \ln(A)}{\partial x} = \kappa_{A} \cosh^{{-}2}{\left( \frac{x-x_0}{w_{A} a}\right)}, \end{equation}

with $T(x_0)=2.14\ \mathrm {keV}$, where $x_0/R_0=0.5$ is the position at which the gradients have their maximum. The density at position $x_0$ is adjusted to obtain the desired value of the reference plasma $\beta (x_0)=8{\rm \pi} n(x_0)T(x_0)/B_{\textrm {ref}}^2$. The parameters $\kappa _{A}$ and $w_{A}$ set the amplitude and the width of the gradients, respectively, and are chosen to be $\kappa _{T_i}=\kappa _{T_e}=6.96$, $w_{T_i}=w_{T_e}=0.3$, $\kappa _{n_i}=\kappa _{n_e}=2.23$ and $w_{n_i}=w_{n_e}=0.3$. The box lengths are chosen to be $(L_{x},L_{v_{||}},L_{\mu })=(80\, \rho _{s},3\,v_{\textrm {th},\sigma } ,9\, T_{0,\sigma }(x_0)/B_{\textrm {ref}})$ and $L_{y}=21.13\, \rho _{s}$, which results in resolving multiples of the toroidal mode number $n_0=19$. Here, $\rho _{s}=c_{s}/\varOmega _{s}$ is the ion sound Larmor radius at reference position $x_0$. The simulation is performed using deuterium as ionic species, as well as electrons that have twice their realistic mass. The resolution at which GENE-3D converged is $(1344\times 16\times 16\times 64\times 32)$ in $(x,y,z,v_{||},\mu )$, respectively, and hyperdiffusion was set as $\eta _{x}=2.0$, $\eta _{y}=0.05$, $\eta _{z}=2.0$ and $\eta _{v_{||}}=0.2$. As is mentioned in Appendix A, periodic boundary conditions are employed in the bi-normal direction, whereas zero-valued Dirichlet boundary conditions are used for the radial domain. In order to avoid numerical instabilities caused by the fixed boundaries, buffer zones with a Krook damping operator are used at the radial boundaries, each covering 5 % of the domain.

Figure 1(a) shows the growth rate $\gamma$ of the dominant linear instability as a function of $\beta$, whereas figure 1(b) shows the frequency $\omega$ of the respective mode. The red lines indicate the results obtained from GENE, whereas blue represents the results of GENE-3D. Excellent agreement is found for both frequencies and growth rates as the relative difference between the results of both codes is below 3 % in all cases. In particular, one can observe the expected stabilisation of the ITG mode by increasing $\beta$ for both codes, as the growth rate decreases for $\beta$ between 0.5 % and 1.4 %, while the frequency changes only slightly in the same interval. For values of $\beta$ between 1.4 % and 1.75 % one can observe a transition of dominant modes from an ITG mode to a kinetic ballooning mode (KBM), which can be identified by the rapid increase in the mode frequency. The results for this interval are only shown here for the GENE code, since resolving this transition requires a much higher resolution than the one that was already used, which makes it impractical to investigate it with GENE-3D. Increasing $\beta$ further results in an amplification of the KBM growth rate, while simultaneously causing a moderate decrease in its frequency, as observed with both codes. Finally, the radial and poloidal structures of the KBM at $\beta =2.5\,\%$ are compared in figure 2, again showing excellent agreement between both codes.

Figure 1. Linear growth rates (a) and mode frequencies (b) of the $n_0=19$ mode as a function of $\beta$. Red shows the results obtained from GENE, blue those from GENE-3D.

Figure 2. Normalised squares of the electrostatic and parallel vector potential for the scenario using $\beta =2.5\,\%$. Radial (a) and poloidal (b) structures of the electrostatic potential, and radial (c) and poloidal (d) structures of the parallel vector potential. The red dashed line shows the results obtained from GENE, the blue solid line those from GENE-3D.

3.2. Nonlinear turbulence at finite plasma-$\beta$ in a tokamak

While the geometry being used is kept the same, it is beneficial for a nonlinear comparison to choose a broader profile than the one used in the previous subsection, as there is less unwanted profile relaxation in a gradient-driven approach with GENE-3D. Therefore, new plasma profiles of the form

(3.4)\begin{equation} \frac{A}{A(x_0)}=\left[ \frac{\cosh\left(\dfrac{x-x_0+\varDelta_{A}}{w_{A}}\right)}{\cosh\left(\dfrac{x-x_0-\varDelta_{A}}{w_{A}}\right)} \right]^{-({\kappa_{A} w_{A}}/{2})({a}/{R_0})} \end{equation}

are chosen, with $x_0/R_0=0.5$, $\kappa _{T_i}=\kappa _{T_e}= 6.66$, $\kappa _n=2.20$, $w_{T_i}=w_{T_e}=w_{n}=0.04$ and $\varDelta _{T_i}=\varDelta _{T_e}=\varDelta _{n}=0.8$. Furthermore, the plasma pressure is such that $\beta _{e}(x_0)=0.75\,\%$. In order to save computational time, the simulations were performed further increasing the electron-to-ion mass ratio to $m_{e}/m_{i}=0.01$ and a finite-size parameter $\rho ^*_{s}=\rho _{s}/L_{\textrm {ref}}=\rho _{s}/R_0= 0.01$. The simulation covers the radial domain $0.1\leqslant x/R_0 \leqslant 0.9$, using buffer zones of 10 % at each side, with normalised box sizes being $(L_{x},L_{y},L_{v_{||}},L_{\mu })=(80.0\,\rho _{s},111.4 \,\rho _{s},3\, v_{\textrm {th},\sigma },9\, T_{0,\sigma }(x_0)/B_{\textrm {ref}})$ and a resolution of $(n_{x},n_{k_y},n_{z},n_{v_{||}},n_{\mu })=(160\times 64\times 24\times 64\times 24)$ for GENE and $(n_{x},n_{y},n_{z},n_{v_{||}},n_{\mu })=(160\times 256\times 24\times 64\times 24)$ for GENE-3D. Finally, heat and particle sources, which are explained in Appendix A, were added with amplitudes being $\kappa _{H}=\kappa _{P}=0.1$ and hyperdiffusion was set as $\eta _{x}=\eta _{y}=0.05$, $\eta _{z}=2.0$ and $\eta _{v_{||}}=0.2$.

In order to compare the results of the two nonlinear simulations, the heat fluxes $Q_{\sigma }$ for species $\sigma$ are considered in the following. The heat flux itself can be split into an electrostatic and an electromagnetic component, $Q_{\sigma }=Q_{\textrm {es},\sigma } +Q_{\textrm {em},\sigma }$, which are defined as

(3.5)\begin{equation} \frac{Q_{\textrm{es},\sigma}}{Q_{\textrm{GB}}} =\frac{m_{\sigma} c}{2 B_0^2} \int v^2\, F_{1,\sigma}^{\textrm{particle}}(\boldsymbol{x},{\boldsymbol{v}},t) \boldsymbol{B}_0\times \boldsymbol{\nabla} \mathcal{G}\left\{\phi_1 \right\} \textrm{d}^3v \end{equation}

and

(3.6)\begin{equation} \frac{Q_{\textrm{em},\sigma}}{Q_{\textrm{GB}}} ={-}\frac{m_{\sigma}}{2c} \int v^2\, F_{1,\sigma}^{\textrm{particle}}(\boldsymbol{x},{\boldsymbol{v}},t) v_{||} \boldsymbol{\nabla} A_{1,||} \times \boldsymbol{b}_0 \, \textrm{d}^3v, \end{equation}

respectively, with

(3.7)\begin{equation} Q_{\textrm{GB}}=n_{\textrm{ref}} \, T_{\textrm{ref}}\, c_{s} (\rho_{s}^*)^2. \end{equation}

Here, $F_{1,\sigma }^{\textrm {particle}}$ is the perturbed distribution function of species $\sigma$, evaluated at particle position.

Figure 3 shows the volume-averaged time traces of the ion and electron heat flux contributions. The fluxes are time-averaged over the interval $t\in [100,345]R_0/c_{s}$, which roughly contains the same number of burst, as follows. The time traces are divided in disjoint intervals each approximately three autocorrelation times long. An average over each one of them is performed, the result of which is then used to compute an ensemble mean and variance of the heat flux. The results are presented in table 1. Both simulations are in good agreement with each other, as the relative difference between their respective ensemble mean values is below 10 % for all components.

Figure 3. Time traces of the volume-averaged electron and ion heat fluxes. Yellow and green lines indicate average over the given time interval for GENE and GENE-3D, respectively. Electrostatic heat flux of ions (a) and electrons (b), and electromagnetic heat flux of ions (c) and electrons (d).

Table 1. Time-averaged heat flux contributions of GENE and GENE-3D.

The radial profiles of the heat flux contributions, averaged over a flux surface and in time, shown in figure 4, recover reasonable agreement between both codes as well. Here GENE and GENE-3D produce similar heat flux profiles, where ion and electron contributions are comparable in their magnitude. For both particle species, the fluxes are dominated by the electrostatic contribution, peaking around $x/R_0=0.6$. At this position, GENE calculates heat fluxes of $Q_{\textrm {es},\textrm {ions}}=103.24\, Q_{\textrm {GB}}$ and $Q_{\textrm {es},\textrm {electrons}}=80.09\, Q_{\textrm {GB}}$, whereas GENE-3D gives $Q_{\textrm {es},\textrm {ions}}=104.72\, Q_{\textrm {GB}}$ and $Q_{\textrm {es},\textrm {electrons}}=75.41\, Q_{\textrm {GB}}$. Both results differ by less than $6\,\%$, giving additional confidence in the numerical implementation of GENE-3D. Finally, the spectra of the electrostatic heat fluxes at $x/R_0=0.6$ are compared in figure 5. There are some small deviations between the heat flux spectra of both codes at the smallest scale, which is the result of the different representations of the bi-normal coordinate. Nevertheless, there is excellent agreement of the spectra in the wavenumber interval where most of the transport is located. Therefore, overall the linear as well as nonlinear verification studies of GENE-3D can be considered successful.

Figure 4. Radial profiles of heat flux contributions. Electrostatic heat flux of ions (a) and electrons (b), and electromagnetic heat flux of ions (c) and electrons (d).

Figure 5. The $k_{y}$ spectra of the electrostatic heat fluxes, evaluated at $x/R_0=0.6$.