Verification of a Fully Implicit Particle-in-Cell Method for the $v_{\parallel}$ Formalism of Electromagnetic Gyrokinetics in the XGC Code

In our scheme, electrons are subcycled [12, 59, 44, 60], meaning they are advanced using several small time steps over the interval $n\Delta t\leq t\leq(n+1)\Delta t$. Hence, we need to define the perturbed electric and magnetic fields over the continuous time interval between steps $n$ and $n+1$. Following [44, 45, 46, 47, 48], we take the electric field to be constant in time over the subcycling interval as

where $\tilde{A}_{\parallel}^{n}$ is defined recursively in time by

and is initialized by solving Eq.(1) at timestep $0$. Consistent with the time derivative of the parallel vector potential being constant within the interval, we take $\delta{\bf B}$ to vary linearly in time as [46, 47, 48]:

The gyrocenter positions of marker particles are described using a cylindrical coordinate system $(R,Z,\varphi)$, where $R$ is the major radius, $Z$ is the distance along the cylindrical axis, and $\varphi$ is the toroidal angle. Together with $v_{\parallel}$ and the $\delta f$ particle weight $w$, there are five evolving variables for each marker particle. These can be written in vector form as $\boldsymbol{\zeta}=[R,Z,\varphi,v_{\parallel},w]^{T}$, and marker particle evolution can then be expressed as

where

In this expression, we have defined the vector quantity ${\bf F}=\langle\delta{\bf E}\rangle-\frac{\mu}{q_{s}}\nabla B$ and have assumed spatial variations in $f_{0s}$ to be along the poloidal flux coordinate $\Psi$. Given the discrete-time potentials $\phi^{n}$, $\phi^{n+1}$, $\tilde{A}_{\parallel}^{n}$, and $A_{\parallel}^{n+1/2}$ defined on the mesh, Eq.(13) and Eq.(15) together with the particle interpolation methods allow for the evaluation of the right hand side vector ${\bf G}(\boldsymbol{\zeta},t)$ at all particle locations $\boldsymbol{\zeta}$ within the mesh and all times within the interval $n\Delta t\leq t\leq(n+1)\Delta t$. This allows considerable freedom in choosing an ODE integration method to advance Eq.(16) in time from step $n$ to $n+1$. Here, we choose a standard fourth order Runge-Kutta method (RK4) for both electrons and ions. For subcycled electrons, we divide the interval into an integer number of sub-intervals and apply RK4 successively over the sub-intervals to advance from step $n$ to $n+1$. Adaptive integration schemes such as the one considered in Ref. [44] are well-suited for subcycling and will be explored in the future. Ions are not subcycled and use field values at the center and ends of the time interval to compute the RK4 stages.

By choosing an implicit time discretization, the particles and fields are interdependent at each timestep. Particles are pushed using fields that depend on $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$, yet these potentials are determined from the particle states over the interval $n\Delta t\leq t\leq(n+1)\Delta t$. In particular, the right hand side of Ampère’s law is computed from the time averaged current depositions over the subcycled timesteps between $n$ and $n+1$, and the right hand side of the gyrokinetic Poisson equation is computed from the number density deposition at $n+1$. We require a self-consistent state between the particles and fields at each timestep, which can be expressed in terms of low-dimensional (compared to the degrees of freedom in the particle system) residuals by using the spatially discretized versions of Eqs.(1)–(2). We define the residuals by

For given $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$, evaluation of the residuals involves: pushing the marker particles from step $n$ to $n+1$, as described in the previous subsection; depositing $\langle\delta j_{\parallel s}\rangle^{n+1/2}$ and $\langle\delta n_{s}\rangle^{n+1}$ to the mesh; and evaluating the elliptic operators in Ampère’s law and the gyrokinetic Poisson equation. Self-consistency is expressed as

which represents a nonlinear system of equations for $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$. We note that both particle and field quantities are coupled in Eq.(19). The moments $\langle\delta j_{\parallel s}\rangle^{n+1/2}$ and $\langle\delta n_{s}\rangle^{n+1}$ are implicitly functions of $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$, where the dependence comes in through the particle equations of motion [44]. Furthermore, given $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$ that satisfy Eq.(19), the consistent particle states follow directly by pushing the particles with the resulting fields.

Our iterative solution method involves wrapping Anderson mixing [50] around a preconditioned Picard iteration scheme. Here, we give a brief overview of the preconditioned Picard iteration scheme. We define a state vector for the potentials ${\bf x}=[\phi^{n+1},A_{\parallel}^{n+1/2}]^{T}$ and a residual vector ${\bf r}=[R_{\phi},R_{A}]^{T}$. In this notation, Eq.(19) can be stated as ${\bf r}({\bf x})=0$, and the preconditioned Picard iteration scheme can be written as

where ${\cal P}$ is an invertible operator whose inverse maps the residual vector to a correction vector $\Delta{\bf x}$ and $k$ is the iteration index. We refer to the operator ${\cal P}$ as the preconditioner and note that the consistency of the iterative method does not depend on the particular form of ${\cal P}$. In other words, any invertible preconditioner used in this scheme will produce a correct solution such that ${\bf r}({\bf x})=0$ provided the scheme converges. The choice of preconditioner will, however, determine if (and at what rate) the scheme converges. As mentioned in the previous subsection, each evaluation of ${\bf r}$ on the right hand side of Eq.(20) involves a particle push and deposition to compute updated velocity moments. One preconditioned Picard iteration from $k$ to $k+1$ therefore follows these steps

1. 1.

   Push particles using the fields computed from the potentials in ${\bf x}^{k}$.

2. 2.

   Deposit $\langle\delta j_{\parallel s}\rangle$ and $\langle\delta n_{s}\rangle$ to compute the residuals in ${\bf r}^{k}$.

3. 3.

   Solve ${\cal P}\Delta{\bf x}^{k}={\bf r}({\bf x}^{k})$ to get the corrections to the potentials $\Delta{\bf x}^{k}$.

4. 4.

   Update the potentials by ${\bf x}^{k+1}={\bf x}^{k}+\Delta{\bf x}^{k}$.

Similarly to what has been reported elsewhere for 6D electromagnetic PIC [46, 47, 48], we choose a preconditioner based on a simplified set of implicitly discretized fluid equations that are linearized with respect to changes in $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$. Our starting point is to consider the continuity and momentum equations for electrons, keeping only the terms relevant to dynamics parallel to the background magnetic field. This greatly simplifies the equations we will need to solve when applying ${\cal P}^{-1}$ to the residuals, yet still captures the terms responsible for the fastest timescales in the system. Since ions evolve on a much slower timescale than electrons, we neglect them entirely in the preconditioner. The starting electron fluid equations are

Next, we consider an implicit discretization of these equations which defines quantities at discrete time locations consistent with the quantities in the PIC system. We have

where the right hand sides of these equations can be determined from information known at the previous timestep. Finally, if we consider the linear responses of the fluid moments due to small perturbations in $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$, we have :

In our notation, $\delta$ refers to perturbations from background quantities in the system, whereas $\Delta$ refers to perturbations due to small changes in the discrete-time potentials $\phi^{n+1}$ and $A_{\parallel}^{n+1/2}$. A 4 $\times$ 4 block matrix ${\mathbb{A}}$ can be written using Eq.(25) together with Ampère’s law and the gyrokinetic Poisson equation, which describes the linear response of the coupled moment-potential system. We write

where the operators in each block are $N\times N$ matrices with $N$ the number of mesh vertices, ${\mathbb{I}}$ is the $N\times N$ identity matrix, and ${\mathbb{M}}$ is a mass matrix for the finite element discretization used in XGC over the poloidal planes. We recall that in step 4 from the previous subsection only the potentials are updated at the end of each iteration. Equation (26) represents an augmented system with equations involving the additional unknowns $\Delta n_{e}$ and $\Delta j_{\parallel e}$. These additional unknowns can provide an intuitive way of expressing the system, where the top two rows of ${\mathbb{A}}$ model the responses of the potentials to changes in the moments, and the bottom two rows model the responses of the moments to changes in the potentials. However, the actual solutions obtained for $\Delta n_{e}$ and $\Delta j_{\parallel e}$ do not play a role in the update.

We use the block matrix system to solve for $\Delta{\bf x}^{k}$ in Eq.(20) in two steps:

1. 1.

   Solve ${\mathbb{A}}{\bf y}={\mathbb{W}}{\bf r}^{k}$, where ${\bf y}=[\Delta\phi,\Delta A_{\parallel},\Delta n_{e},\Delta j_{\parallel e}]^{T}$ and

   $${\mathbb{W}}=\left[\begin{array}[]{cc}{\mathbb{I}}&0\\ 0&{\mathbb{I}}\\ 0&0\\ 0&0\end{array}\right].$$

2. 2.

   Restrict the solution to the potential correction vector $\Delta{\bf x}^{k}$ with $\Delta{\bf x}^{k}={\mathbb{W}}^{T}{\bf y}$.

The block matrix is therefore related to the preconditioner operator in Eq.(20) by ${\cal P}^{-1}={\mathbb{W}}^{T}{\mathbb{A}}^{-1}{\mathbb{W}}$. In our implementation, we use the suite of solvers available in PETSc [61, 62, 63] to invert ${\mathbb{A}}$. Some initial convergence results for the iterative scheme applied to the preconditioned implicit PIC equations are given in Sec. 5. However, detailed studies of convergence rates across different parameter regimes, modifications to the fluid-based preconditioner for improving the convergence rates, and methods for efficiently inverting ${\mathbb{A}}$ will be the subject of a separate paper.

For the toroidal system, we again consider the cylindrical coordinate system $(R,Z,\varphi)$ from Sec. 3.1. In addition, we can define a simple toroidal coordinate system by expressing $R$ and $Z$ in terms of the $R$ coordinate at the magnetic axis $R_{0}$, a minor radius coordinate $r$, and a poloidal angle coordinate $\theta$ as

The background magnetic field model used is from [64] and can be expressed in the cylindrical coordinates as

where $\Psi$ is a poloidal flux function. An analytical form is taken for $\Psi$ depending on $r$ only as

where $q(r)$ is the safety factor. Hence, the magnetic equilibrium is determined once we specify the parameters $B_{0}$ and $R_{0}$ and the function $q(r)$. We take $q(r)$ to be a quadratic polynomial in $r/a$, where $a$ is the minor radius

Hence $q(r)$ is determined by the parameters $q_{0}$, $q_{1}$, and $q_{2}$ in addition to the minor radius $a$. To adapt this model to the periodic cylindrical system used in the SAW benchmark, Eq.(29) is modified by replacing $R$ in the denominator of both terms with the constant $R_{0}$. The toroidal angle $\varphi$ then serves as an axial coordinate and $(R,Z)$ are coordinates within planes perpendicular to the cylindrical axis. The distance along the cylindrical axis is parameterized by $R_{0}\varphi$, and the system remains $2\pi$-periodic in $\varphi$. In this way, we can interpret the cylindrical system as a “straightened” toroidal system.

To describe density and temperature profiles, we first define a normalized logarithmic gradient. Letting $A$ represent either density or temperature, the normalized logarithmic gradient of $A$ is defined as

and the profile of $A$ is then

where $r_{0}$ is a reference value of $r$. We choose an analytical form for $\kappa_{A}$ as

Hence profiles are specified by the parameters $r_{0}$, $R_{0}$, $a$, $w_{A}$, $\kappa_{A}(r_{0})$, and $A(r_{0})$.

We simulate SAW propagation in a periodic cylindrical system with uniform temperature, density, and safety factor profiles. For simplicity, ions are modeled by a uniform background density $n_{0}$ in addition to the ion polarization density in the gyrokinetic Poisson equation Eq. (2). Electrons are kinetic. Considering only the polarization response for the ions is sufficient for providing a minimal model supporting Alfvén wave propagation, and similar models have been considered in the past, for example in Refs. [65, 66, 67]. Physically, the ion polarization response provides a cross-field current to counter the parallel electron current in a bid to maintain quasi-neutrality. In this benchmark, deuterium ions are considered along with electrons at twice the realistic mass, giving a mass ratio of $\frac{m_{i}}{m_{e}}=1836$. The relevant fixed parameters for this test case are given in Table 1.

In Fig.(1), we scan over $\beta_{e}=\frac{\mu_{0}n_{0}T_{e}}{B^{2}}$, which we vary by choosing different values of density ranging from $6.25\times 10^{17}\ \mathrm{m}^{-3}$ to $3.0\times 10^{19}\ \mathrm{m}^{-3}$. Simulations are initialized with a sinusoidal electron density perturbation of the form:

where $A$ is a radial envelope function and c.c. denotes the complex conjugate. We note that $\beta_{e}$ scans of SAW propagation have long been used to verify electromagnetic gyrokinetic PIC codes. For example, [68] and [10] consider similar tests in shearless slab geometry. To shed light on potential cancellation issues, we consider a long wavelength large $\beta$ regime. We initialize our simulations taking $(n,m)=(2,2)$. Furthermore, a Fourier filter is used to keep only the $n=2$ mode. A timestep size of $\Delta t=5.6\times 10^{-3}R_{0}/c_{s}$ is used, where $c_{s}$ is the sound speed given by $c_{s}=\sqrt{T_{e}/m_{i}}=3.1\times 10^{5}\ \mathrm{m/s}$.

The mesh in XGC consists of a number of identical $R-Z$ planes equally spaced in $\varphi$, where within each $R-Z$ plane, an unstructured triangular mesh is employed. For this case, 9038 mesh vertices are used within each $R-Z$ plane, corresponding to a spatial resolution of $\Delta R\approx\Delta Z\approx 9.3\times 10^{-3}\ \mathrm{m}$, which is roughly $0.3$ times the size of the ion-sound gyroradius $\rho_{s}=c_{s}/\Omega_{i}=2.8\times 10^{-2}\ \mathrm{m}$. Between the $R-Z$ planes, a resolution of $R_{0}\Delta\varphi\approx 19.5\ \mathrm{m}$ is used. Finally, we take approximately 50 particles per mesh vertex. We compare the measured real frequencies from the simulations to a simple analytic dispersion relation. In B, we present a dispersion analysis for the SAW resulting in

We note that several simplifying assumptions have been made in the derivation of Eq.(36). For example, we have neglected kinetic effects, effects due to nonuniformities in ${\bf B}$, and finite $k_{\perp}$ effects. For $r/R_{0}\ll 1$, with the magnetic field described in Eq.(29) and Eq.(30) and perturbations of the form Eq.(35), we have $k_{\parallel}\approx(n-m/q)/R_{0}$. Figure (1) shows excellent agreement when comparing simulation results to this approximate dispersion relation, demonstrating the absence of cancellation errors in our implementation.

In addition to comparing real frequencies, we also show the time history of the real and imaginary parts of the $(n,m)=(2,2)$ mode amplitude of $A_{\parallel}$ and $\phi$ for the simulation at $\beta_{e}=3.85\%$ in Fig.(2). A $90^{\circ}$ phase shift is observed between $\mathrm{Re}(\phi)$ and $\mathrm{Im}(A_{\parallel})$ and between $\mathrm{Im}(\phi)$ and $\mathrm{Re}(A_{\parallel})$. This phase relation follows from Eq.(1), Eq.(2), and the electron continuity equation as is shown in A.

Next, we consider a cyclone-like [69] test case in toroidal geometry similar to the case considered in [70], demonstrating the ITG-KBM transition. ITG represents ion temperature-gradient driven modes, and KBM represents kinetic ballooning modes. In this study, we compare results from our fully implicit electromagenetic algorithm in XGC to results obtained from GEM, GENE, and explicit XGC. The GEM and GENE codes are both global $\delta f$ gyrokinetic codes using field-line-following coordinates. GEM is a particle-in-cell code, and GENE is Eulerian using a spectral method in the binormal direction. The electromagnetic capability in GEM is based on the $p_{\parallel}$-formalism, and a modified mass matrix is used in Ampère’s law to mitigate accuracy problems at high $\beta$ from the cancellation problem [11]. In GENE, $v_{\parallel}$ is used as a coordinate for electromagnetic simulations, and a transformation of the perturbed distribution function involving $A_{\parallel}$ is performed to eliminate the time derivative in the inductive component of the electric field [22]. GEM and GENE both use forms of the polarization density in the gyrokinetic Poisson equation that are valid for arbitrary wavelengths. In this study, on the other hand, the implicit and explicit electromagnetic XGC have the short wavelength Padé correction term turned off, as shown in Eq.(2). This may cause some discrepancy for short wavelength modes, which has been shown, however, to be insignificant in the recent $n=19$ benchmarking exercise between GENE and explicit electromagnetic XGC in the Görler benchmarking plasma [41]. Our benchmarking is at the toroidal mode number $n=6$. More detailed descriptions of GEM and GENE can be found in the literature including Refs. [23, 22, 71, 21, 72] for GENE and Refs. [9, 10, 11] for GEM. For further comparison, we mention that XGC uses unstructured triangular meshes in cylindrical coordinates. In addition, the typical mode of operation in XGC is the total-$f$ method described in [6, 7]. However, in this paper XGC is run in standard $\delta f$ mode for a proper comparison to GEM and GENE.

In this benchmark problem, we take a reduced magnetic field strength compared to the case considered in [70]. The motivation for this choice is to reduce the resolution requirements, allowing this benchmark problem to be run at a lower computational cost. In [70], the size parameter $\rho^{*}$, defined as the ratio of ion gyroradius to the minor radius, was approximately $1/180$; here, $\rho^{*}\approx 1/50$. Since the perpendicular resolution requirement is set by the size of the ion gyroradius, the number of mesh nodes required for a well-resolved simulation of this benchmark problem is roughly a factor of $(180/50)^{2}\approx 13.0$ smaller than that required for the case considered in [70]. Both ions and electrons are treated kinetically. We take hydrogen ions and electrons with realistic mass, giving a mass ratio of $\frac{m_{i}}{m_{e}}=1836$, and we take $T\equiv T_{i}=T_{e}$. Furthermore, only the $n=6$ mode is kept in the simulations. We note that the ITG-KBM transition case in [70] used $n=19$. Fixed parameters for defining the magnetic geometry and profiles are given in Table 2.

In Fig.(3), we show the density and temperature profiles normalized by their values at $r_{0}$ on the left and the normalized logarithmic gradients of density and temperature on the right. Compared to the profiles given in Figure 2 of [70], the values of $\kappa_{T}$ and $\kappa_{n}$ are nearly identical at $r_{0}$; however, the $\kappa$ profiles considered here are much flatter in the center and nearly zero at the ends. This results in flat density and temperature profiles near the magnetic axis and the outer radial boundary. More localized gradients are used in this benchmark problem since this can help reduce effects at the boundaries of the simulation domain, which can be more significant in problems with larger $\rho^{*}$.

In Fig.(4), the safety factor and magnetic shear profiles are shown. Note that these are identical to the ones used in [70].

We again scan over $\beta_{e}$ by varying the density value of density at $r_{0}$ from $6.5\times 10^{16}\ \mathrm{m}^{-3}$ to $1.625\times 10^{18}\ \mathrm{m}^{-3}$. In Fig.(5), we compare measured real frequencies and growth rates obtained from the implicit version of XGC to those obtained in GEM, GENE, and explicit XGC. The real frequencies are given in table 3 and the growth rates in table 4 for the cases that were simulated by all four codes/algorithms. Good agreement is observed with each code/algorithm finding the critical value of $\beta_{e}$ to be between $0.65\%$ and $0.75\%$. Given the significant differences that exist between the algorithms and formulations used, we consider the overall agreement in real frequencies and growth rates for this benchmark case to be very good. We note that all four comparisons found a collisionless trapped electron mode (CTEM) to be present at $\beta_{e}=0.65\%$, characterized by a change in the direction of propagation from the ion diamagnetic direction to the electron diamagnetic direction. We note that there is some disagreement between GENE and the other codes in the real frequency for the CTEM case. A possible cause for this disagreement may be the form of the curvature drift that was used in GENE for this benchmarking exercise. GENE was run using the form of the curvature drift given in C, which is different than what was used in XGC and GEM. Both versions of XGC and GEM used the ${\bf B}^{*}$ term in Eq.(5), which implicitly contains the curvature drift. When a true MHD equilibrium is considered, both forms of the curvature drift should be equivalent to first order in the gyrokinetic ordering. However, MHD equilibrium is not ensured in our simplified analytic models of the magnetic field and profiles as in the Görler benchmarking case [70]. Since the physics of the CTEM mode depends strongly on the magnetic curvature, it may be more sensitive to differences in this term than the ITG or KBM modes.

Finally, we show the developed mode structures for $\phi$ and $A_{\parallel}$ in Fig.(6) for the $\beta_{e}=0.05\%$ ITG case and in Fig.(7) for the $\beta_{e}=1.25\%$ KBM case. There is good overall qualitative agreement in the mode structures between the various codes/algorithms for both cases. We note that the electrostatic potentials produced from GEM and GENE feature fine-scale radial structures near rational surfaces. The two algorithms in XGC, on the other hand, produce smoother electrostatic potentials in the radial direction. This may be due to the long-wavelength form of the ion polarization density used in Eq.(2) for the XGC algorithms.

The numerical resolutions used by each code/algorithm to obtain these results are the following. For the fully implicit version of XGC, a timestep size of $\Delta t=4.0\times 10^{-2}\ \frac{L_{T_{i}}}{c_{s}}$ was used, where $L_{T_{i}}=\frac{R_{0}}{\kappa_{T_{i}}(r_{0})}$. Here, 28439 mesh vertices were used in the $R\!\!-\!\!Z$ planes, corresponding to a spatial resolution of approximately $0.5\ \rho_{i}$. Along the toroidal direction, 8 $R\!\!-\!\!Z$ planes were taken within a wedge spanning $1/6$ of a full torus, and 50 particles per mesh vertex are used. The explicit version of XGC used the same mesh and toroidal resolution as was used in the implicit version. The timestep size was $\Delta t=1.6\times 10^{-3}\ \frac{L_{T_{i}}}{c_{s}}$ and 25 particles per mesh vertex were used. GEM used a radial resolution of $0.31\ \rho_{i}$ and resolution in the binormal direction of $0.29\ \rho_{i}$. Along the field line, 64 grid points were used spanning $-\pi<\theta\leq\pi$. The timestep size used in GEM was $\Delta t=5.0\times 10^{-3}\ \frac{L_{T_{i}}}{c_{s}}$ and 16 marker ions and 32 marker electrons per cell were used. Finally, GENE used a radial resolution of $0.19\ \rho_{i}$ and used a single mode in the binormal direction corresponding to a toroidal mode number of $n=6$. Along the field line, 32 grid points were used spanning $-\pi<\theta\leq\pi$. In velocity space, 48 grid points were used for $v_{\parallel}$ and 24 grid points were used for $\mu$. The timestep size used in GENE was $\Delta t=8.8\times 10^{-3}\ \frac{L_{T_{i}}}{c_{s}}$.

Each code used homogeneous Dirichlet boundary conditions when solving for $\phi$ and $A_{\parallel}$. However, there was some variation in the locations of the radial boundaries. The implicit version of XGC took only an outer radial boundary at $r/a=0.9$. Since XGC uses cylindrical coordinates, it is able to include the magnetic axis within the simulation domain. The explicit version of XGC took an outer radial boundary at $r/a=0.98$ and an inner radial boundary at $r/a=0.07$. In addition, filtering was applied after the field solves to set the potentials outside of the region $0.24\leq r/a\leq 0.83$ to zero. GEM took an outer radial boundary at $r/a=0.9$ and an inner radial boundary at $r/a=0.1$. GENE took an outer radial boundary at $r/a=0.975$ and an inner radial boundary at $r/a=0.025$. Both versions of XGC and GEM set particle weights to zero when particles exited the simulation domain, and GENE used a buffer region outside the simulation domain with a Krook collision operator designed to damp the perturbed distribution function towards zero.

For completeness, we include here some timing data from each code/algorithm in simulating the $\beta_{e}=1.25\%$ case. Since the runs using the explicit version of XGC in the benchmarking study were carried out without consideration of performance, we use a more practical timestep size here for a fairer comparison. All simulations were run using 16 Haswell processor nodes on the Cori supercomputer at NESRC. Both versions of XGC and GENE ran using MPI parallelization with a total of 1024 MPI processes. GEM ran using a hybrid MPI/OpenMP parallelization with 128 MPI processes and 8 OpenMP threads per MPI process. We choose our metric to be the wallclock time required to simulate a physical time unit of $L_{T_{i}}/c_{s}$. GENE required a wallclock time of 2.4 seconds, GEM required 75.0 seconds, the implicit version of XGC required 110.3 seconds, and the explicit version of XGC required 176.6 seconds. We note that in collecting this specific timing data, the explicit version of XGC used a more practical timestep size of $1.0\times 10^{-2}\frac{L_{T_{i}}}{c_{s}}$ and 50 particles per mesh vertex, in an attempt to improve the unnecessarily small timestep size used in the physics benchmarking exercise in the explicit version of XGC and to make the particle number the same. All other parameters are the same as those used in the benchmarking exercise presented above in this section. We caution here that the timing data comes from a single common simulation scenario in which parameters were chosen to provide a low cost physics benchmarking problem without considering computational algorithm optimization, number of subcycling steps (4 steps used in this study is at the low end), optimal hardware conditions, programing language, and practical problem size per each code. For example, optimization of the production XGC has been known to scale well for big physics problems on extreme scale computers that are equipped with high-performance accelerators. In addition, optimal choices of parameters for balancing both accuacy and efficiency are yet to be determined for the different versions of XGC. Hence, the performance numbers provided here should not be used to draw any conclusions about the efficiency of the codes/algorithms in simulating different realistic problems targeted by each codes.