# Kinetic Simulation of Collisional Magnetized Plasmas with Semi-Implicit Time Integration

Plasmas with varying collisionalities occur in many applications, such as tokamak edge regions, where the flows are characterized by significant variations in density and temperature. While a kinetic model is necessary for weakly-collisional high-temperature plasmas, high collisionality in colder regions render the equations numerically stiff due to disparate time scales. In this paper, we propose an implicit-explicit algorithm for such cases, where the collisional term is integrated implicitly in time, while the advective term is integrated explicitly in time, thus allowing time step sizes that are comparable to the advective time scales. This partitioning results in a more efficient algorithm than those using explicit time integrators, where the time step sizes are constrained by the stiff collisional time scales. We implement semi-implicit additive Runge-Kutta methods in COGENT, a finite-volume gyrokinetic code for mapped, multiblock grids and test the accuracy, convergence, and computational cost of these semi-implicit methods for test cases with highly-collisional plasmas.

## Authors

• 2 publications
• 2 publications
• 2 publications
• 2 publications
• ### Convergence estimates of a semi-Lagrangian scheme for the ellipsoidal BGK model for polyatomic molecules

In this paper, we propose a new semi-Lagrangian scheme for the polyatomi...
02/29/2020 ∙ by Sebastiano Boscarino, et al. ∙ 0

• ### Parallel-in-Time Solution of Eddy Current Problems Using Implicit and Explicit Time-stepping Methods

The time domain analysis of eddy current problems often requires the sim...
12/14/2020 ∙ by Idoia Cortes Garcia, et al. ∙ 0

• ### An all speed second order well-balanced IMEX relaxation scheme for the Euler equations with gravity

We present an implicit-explicit well-balanced finite volume scheme for t...
12/19/2019 ∙ by Andrea Thomann, et al. ∙ 0

• ### IMEX HDG-DG: a coupled implicit hybridized discontinuous Galerkin (HDG) and explicit discontinuous Galerkin (DG) approach for shallow water systems

We propose IMEX HDG-DG schemes for planar and spherical shallow water sy...
11/07/2017 ∙ by Shinhoo Kang, et al. ∙ 0

• ### Exponential Time Differencing for the Tracer Equations Appearing in Primitive Equation Ocean Models

The tracer equations are part of the primitive equations used in ocean m...
10/05/2019 ∙ by Sara Calandrini, et al. ∙ 0

• ### Hierarchical Optimization Time Integration for CFL-rate MPM Stepping

We propose Hierarchical Optimization Time Integration (HOT) for efficien...
11/18/2019 ∙ by Xinlei Wang, et al. ∙ 0

• ### A stable semi-implicit algorithm

When the singular values of the evolution operator are all smaller or al...
05/11/2019 ∙ by João P. S. Bizarro, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The purpose of this paper is to describe the application and performance of a semi-implicit time integration algorithm for the solution of a system of Vlasov-Fokker-Planck equations, motivated by the goal of simulating the edge plasma region of tokamak fusion reactors. Plasma dynamics in the tokamak edge region is an unsteady multiscale phenomenon, characterized by a large range of spatial and temporal scales due to the density and temperature variations. Figure 1 shows the cross-section of a typical tokamak fusion reactor. The geometry is defined by the magnetic flux surfaces that contain the plasma, and the core and the edge regions are marked. Within the edge region, as the temperature decreases from the hot near-core region to the cold outer-edge region, there are three scale regimes. The hot and dense plasma in the inner edge region adjacent to the core plasma is weakly collisional, and the mean free paths of the particles are significantly larger than the density and temperature gradient length scales cohenxu2008 ; dorf_etal_2013 ; dorf_etal_2014 . Near the separatrix that separates the closed magnetic field lines from the open ones, the plasma is moderately collisional, and the mean free paths are comparable to the density and temperature gradient scales. At the outer edge (near the material surfaces), the cold plasma is strongly collisional. In this region, the particle mean free paths are significantly smaller than the density and temperature length scales.

Due to the weak collisionality near the core, a kinetic description with an appropriate collision model is required to model accurately the perturbations to the velocity distribution from the Maxwellian distribution. The Vlasov-Fokker-Planck (VFP) equation governs the evolution of the distribution function of each charged particle in the position and velocity space thomasetal2012 . In the presence of a strong, externally-applied magnetic field, the ionized particles gyrate around the magnetic field lines. In the context of tokamaks, the radius of this gyromotion (gyroradius) is much smaller than the characteristic length scales. The gyrokinetic VFP equation represents the dynamics of the particles averaged over the gyromotion lee1983 ; cohenxu2008 , i.e., it describes the motion of the guiding center of the particles. It is thus expressed in terms of the gyrocentric coordinates (parallel and perpendicular to the magnetic field lines) and does not contain phase-dependent terms. As a result, one velocity dimension and the fast time scale of the gyromotion is removed. The drift-kinetic model is the long wavelength limit of the gyrokinetic model, where turbulent length scales (comparable to the gyroradius) are neglected. We consider the drift-kinetic VFP equation in this paper.

Numerical algorithms to solve the gyrokinetic equation can be classified into three families. The Lagrangian particle-in-cell (PIC) approach

changku2008 ; heikkinenetal2001 ; heikkinenetal2006 ; chenetal2011 ; chenchacon2014 ; chenchacon2015 solves the equations of motion for “superparticles” in the Lagrangian frame. The primary drawback is that the number of superparticles needed to control the numerical noise is very large thomasetal2012 ; filbetetal2001 . Eulerian methods idomuraetal2008_2 ; scott2006 ; xuetal2007 ; cohenxu2008 ; xuetal2008 ; xuetal2010 solve the VFP equation on a fixed phase space grid, where the phase space comprises the spatial position and velocity coordinates. Semi-Lagrangian approaches gysela ; rossmanithseal2011 ; crouseillesetal2009 ; qiuchristlieb2010 ; qiushu2011

follow characteristics, either tracing backward or forward in time, requiring interpolation of the solution from or to a fixed phase-space grid. In this paper, we adopt the Eulerian approach that allows the use of advanced numerical algorithms developed by the computational fluid dynamics (CFD) community.

The gyrokinetic VFP equation is a parabolic partial differential equation (PDE) that comprises an advective Vlasov term and the Fokker-Planck collision term. The Vlasov term describes incompressible flow in the phase space

thomasetal2012 ; hahm1996 , i.e., the evolution of the distribution function by the gyrokinetic velocity of the particles and their acceleration parallel to the magnetic field. The Fokker-Planck term describes the collisional relaxation of the distribution function to the Maxwellian distribution thomasetal2012 ; rosenbluth1957 and acts in the velocity space only. The difference in the time scales of these two terms varies significantly in the edge region due to the variation in the collisionality. The weakly and moderately collisional plasma at the inner edge and the separatrix, respectively, exhibits collisional time scales much smaller or comparable to the Vlasov time scales, while the strongly-collisional plasma at the outer edge exhibits collisional time scales much smaller than the Vlasov time scales. Consequently, the VFP equation exhibits disparate scales and is numerically stiff. Explicit time integration methods are inefficient because the time step size is constrained by the collisional time scale. Implicit time integration methods, on the other hand, are unconditionally linearly stable, but they require solutions to a linear or nonlinear system of equations. Linearly implicit lemou2005 ; buetetal1997 ; lemou1998 ; larroche1993 ; larroche2003 ; lemou2004 ; casanovaetal1991 and nonlinearly implicit methods larsenetal1985 ; epperlein1994 ; chaconetal2000_2 ; taitanoetal2015 ; taitanoetal2016 have been applied to the Fokker-Planck operator in isolation, and several implicit algorithms have also been proposed for the VFP equation epperleinetal1988 ; mousseauknoll1997 ; kinghambell2004 ; thomasetal2009 . In our context, the accurate resolution of the unsteady dynamics requires time steps comparable to the Vlasov time scales, and thus a fully-implicit approach is inefficient.

Semi-implicit or implicit-explicit (IMEX) time integration methods ascherruuthspiteri ; pareschirusso ; kennedycarpenter allow the partitioning of the right-hand-side (RHS) into two parts: the stiff component integrated implicitly, comprising time scales faster than those of interest, and the nonstiff component integrated explicitly, comprising the slower scales. Such methods have been applied successfully to many multiscale applications durranblossey ; giraldorestellilauter2010 ; ghoshconstaSISC2016 ; jacobshesthaven2009 ; kumarmishra2012 ; kadiogluetal2010 ; salmietal2014 ; wolkeknoth2000 . In the context of the Vlasov equation, semi-implicit approaches have been proposed chengknorr1976 ; candywaltz2003 ; idomuraetal2008_1 ; maeyamaetal2012 that partition the Vlasov term and integrate the parallel advection implicitly in time.

In this paper, we propose a semi-implicit algorithm for the gyrokinetic VFP equation, where we integrate the Vlasov term explicitly while the stiff Fokker-Planck term is integrated implicitly to allow time steps that are comparable to the Vlasov time scales. The Fokker-Planck term represents a nonlinear advection-diffusion operator that can be expressed in two forms: the Landau form landau1937 and the Rosenbluth form rosenbluth1957 . Although the two forms are theoretically equivalent, they differ in their implications on the implicit numerical solution of this term. The Landau form expresses the advection and diffusion coefficients as direct integrals and is well-suited for conservative numerical methods chang1970 ; kho1985 ; larsenetal1985 ; epperlein1994 ; lemou2005 ; buetcordier1998 ; berezinetal1987 . However, when integrated implicitly, the integral form results in a dense, nonlinear system of equations. Naive algorithms scale as , where is the number of velocity space grid points, whereas fast algorithms filbet2002 ; buetetal1997 ; banksetal2016 ; brunneretal2010 scale as . The Rosenbluth form, on the other hand, relates the advection and diffusion coefficients to the distribution function through Poisson equations for the Rosenbluth potentials in the velocity space. Thus, each evaluation of the Fokker-Planck term requires the solution to the Poisson equations, which scales as , assuming that an efficient Poisson solver is available. The primary difficulty is the need for solving the Poisson equations, defined on an infinite velocity space, on the truncated numerical domain, and several approaches have been proposed james1977 ; bellicandy2012 ; mccoyetal1981 ; landremanernst2013 ; pataki2011 ; larroche1993 ; chaconetal2000_2 . An additional difficulty with the Rosenbluth form is the difficulty in enforcing mass, momentum, and energy conservation chaconetal2000_1 ; taitanoetal2015 , and modifications are necessary to ensure energy and momentum are conserved to round-off errors. In this paper, we consider the Rosenbluth form of the Fokker-Planck term.

Our semi-implicit approach is implemented in COGENT dorf_etal_2012 ; dorf_etal_2013 ; dorf_etal_2014 , a high-order finite-volume code that solves the gyrokinetic VFP equations on mapped, multiblock grids representing complex geometries. We implement multistage, conservative additive Runge-Kutta (ARK) methods kennedycarpenter , where the resulting nonlinear system of equations for the implicit stages are solved using the Jacobian-free Newton-Krylov approach knollkeyes2004 . We investigate the performance of the ARK methods for test problems representative of the tokamak edge and compare it to that of the explicit Runge-Kutta (RK) methods for VFP problems. In particular, we verify the accuracy and convergence of the ARK methods and their computational efficiency with respect to the RK methods. Although the preconditioning that we have implemented reduces the cost of the implicit solve significantly, our approach is a preliminary effort, and the implementation of better preconditioning techniques will be addressed in the future.

The outline of the paper is as follows. Section 2 introduces the gyrokinetic VFP equations and the nondimensionalization used in our implementation. Section 3 describes the high-order finite-volume method that COGENT uses to discretize the equations in space. The time integration methods are discussed in Section 4. The algorithm verification through two test cases is presented in Section 5, and Section 6 summarizes the contributions of this paper.

## 2 Governing Equations

The plasma dynamics at the tokamak edge are described by the full- gyrokinetic VFP equation hahm1996 ; dorf_etal_2013 . In this paper, we consider a single species, and the axisymmetric governing equation is expressed in the nondimensional form as

 ∂B∗∥f∂t+∇R⋅(˙RB∗∥f)+∂∂v∥(˙v∥B∗∥f)=cf(f), (1)

where is the distribution function defined on the phase space .

is the spatial gyrocenter position vector in the configuration space with

as the radial coordinate and as the poloidal coordinate, is the velocity parallel to the externally applied magnetic field , and

is the magnetic moment. The configuration space

is two-dimensional, and (1) is a four-dimensional (2D–2V) PDE. denotes the divergence operator in the configuration space. In the long wavelength limit, the gyrokinetic model reduces to the drift-kinetic model, and the Vlasov velocity and parallel acceleration are given by

 ˙R ≡˙R(R,v∥,μ,t)=1B∗∥[v∥B∗+λaZ^b×(ZE+μ2∇RB)], (2a) ˙v∥ ≡˙v∥(R,v∥,μ,t)=−1mB∗∥B∗⋅(ZE+μ2∇RB), (2b)

where and are the mass and ionization state, respectively, is the Larmor number (ratio of the gyroradius to the characteristic length scale), and are the unit vector along the magnetic field and magnitude of the magnetic field, and

 B∗≡B∗(R,v∥)=B+λamv∥Z∇×^b. (3)

is the Jacobian of the transformation from the lab frame to the gyrocentric coordinates. In this paper, we consider cases with uniform magnetic field, and thus . The electric field is , where is the electrostatic potential. In this paper, we consider cases where is specified; however, COGENT generally computes a self-consistent electrostatic potential from the species charge densities dorf_etal_2012 ; dorf_etal_2013 where a Boltzman or vorticity model is assumed for the electrons.

The Rosenbluth form of the Fokker-Planck collision term for a single species is expressed as dorf_etal_2014

 c(f) =νc∇(v∥,μ)⋅→Γ, (4a) →Γ =→σf+↔κ∇(v∥,μ)f, (4b)

where is the collision frequency, and the coefficients are defined as follows:

 (5)

where

 σv∥=∂φ∂v∥,   σμ=4μ(mB)∂φ∂μ, κv∥v∥=−∂2ϱ∂v2∥, κv∥μ=κμv∥=−4μ(mB)∂2ϱ∂v∥∂μ, κμμ=−8μ(mB)2[2μ∂2ϱ∂μ2+∂ϱ∂μ]. (6)

The operator denotes the gradient operator in the velocity space . The Rosenbluth potentials, and , are related to the distribution function through the following Poisson equations in the two-dimensional velocity space :

 ∂2φ∂v2∥+mB∂∂μ(2μ∂φ∂μ) =f, (7a) ∂2ϱ∂v2∥+mB∂∂μ(2μ∂ϱ∂μ) =φ. (7b)

The Maxwellian distribution function, defined as

 fM(R,v∥,μ)=n(R)√π(m2T(R))32exp⎛⎝−m{v∥−¯¯¯¯¯v∥(R)}2+μB(R)2T(R)⎞⎠, (8)

is an exact solution of the Fokker-Planck collision operator (), where the temperature is given by

 T(R) =1n(R)∫∫13[(v∥−¯¯¯¯¯v∥(R))2+μB]f(R,v∥,μ)B∗∥(R,v∥)dv∥dμ (9)

and

 n(R) =1m∫∫f(R,v∥,μ)B∗∥(R,v∥)dv∥dμ, (10) ¯¯¯¯¯v∥(R) =1mn(R)∫∫v∥f(R,v∥,μ)B∗∥(R,v∥)dv∥dμ (11)

are the number density and average parallel velocity, respectively.

The equations in the preceding discussion involve non-dimensional variables. The non-dimensionalization is derived by

 t=~ttref,  v∥=~v∥vref,  n=~nnref,  m=~mmref,  f=~ffref,  T=~TTref, B=~BBref,  ϕ=~ϕϕref,  μ=~μμref,    x,y,L=~x,~y,~LLref,

where is the physical (dimensional) quantity and are the reference quantities. The Larmor number and the collision frequency are expressed as

 λa=mrefvrefeLrefBref,     νc=freftrefe4m2refΛc(4πZ2m)2,

respectively, where is the Coulomb logarithm. The primitive reference parameters are (number density), (temperature), (length), (mass), and (magnetic field); the derived reference quantities are

 vref=√Trefmref (thermal speed),  tref=Lrefvref (transit time), μref=Tref2Bref (magnetic moment),  fref=nrefπv3ref (distribution function), ϕref=Trefe (potential);

and is the elementary charge.

## 3 Spatial Discretization

Equation (1) is a parabolic PDE, which we discretize in space using a fourth-order finite-volume method on a mapped grid collelaetal2011 ; mccorquodaleetal2015 . The computational domain is defined as a four-dimensional hypercube of unit length in each dimension,

 Ω={ξ:0≤ξd≤1,1≤d≤D};  ξd=ξ⋅ed, (12)

partitioned by a uniform grid with a computational cell defined as

 \upomegai=D∏d=1[(i−12ed)h,(i+12ed)h], (13)

where is the number of spatial dimensions, is the unit vector along dimension , and is an integer vector representing a four-dimensional grid index. The transformation defines the mapping between a vector in the physical space and a vector in the computational space . Equation (1) is integrated over a physical cell to yield the integral form

 (14)

where the Vlasov and collision terms are

 V(B∗∥f)≡−[∇R⋅(˙RB∗∥f)+∂∂v∥(˙v∥B∗∥f)],  C(f)=cf(f). (15)

Defining the computational cell-averaged solution as

 ¯fi=1\upomegai–––––––––––∫\displaylimits\upomegaiB∗∥fJdξ;   J≡∣∣∣∂x∂ξ∣∣∣, (16)

where is the volume of the computational cell , the physical cell-averaged solution is

 ˘fi=1X(\upomegai)––––––––––––––––∫\displaylimitsX(\upomegai)B∗∥fdx=(∫\displaylimitsViJdξ)−1∫\displaylimits\upomegaiB∗∥fJdξ=¯J−1i¯fi, where ¯Ji=1\upomegai–––––––––––∫\displaylimits\upomegaiJdξ. (17)

The Vlasov and collision terms in (14) can be written in the divergence form as follows:

 V(B∗∥f) =∇x⋅V(B∗∥f);  V(B∗∥f)=−[(˙R⋅^r)B∗∥f,(˙R⋅^θ)B∗∥f,˙v∥B∗∥f,0]T, (18a) C(f) =∇x⋅C(f);  C(f)=νc[0,0,→Γ⋅^v∥,→Γ⋅^μ]T, (18b)

where are the unit vectors along the coordinates, respectively, and denotes the divergence operator in the physical space. Using the divergence theorem, (14) is expressed as

 ∂¯fi∂t =1\upomegai–––––––––––[∫\displaylimits∂X(\upomegai)V(B∗∥f)⋅^nds+∫\displaylimits∂X(\upomegai)C(f)⋅^nds] =1\upomegai–––––––––––D∑d=1[∫\displaylimitsAi+12edNT(V+C)dAξ−∫\displaylimitsAi−12edNT(V+C)dAξ], (19)

where is the outward normal for , , and denote the faces along dimension of the cell . The face-averaged Vlasov and collision fluxes on a given face are defined as

 (20a) (20b)

The computational cell is a hypercube with length in each dimension with the volume as , and the area of each face as . Thus, (19) can be written as

 (21)

Defining the solution vector as consisting of the cell-averaged distribution function at all the computational cells in the grid, (21

) is expressed for the entire domain as a system of ordinary differential equations (ODEs) in time,

 d¯fdt (22a) ^V(¯fi) =1hD∑d=1(⟨^Vi+12ed⟩−⟨^Vi−12ed⟩)=V(B∗∥f)+O(ΔRp,Δvp∥,Δμp), (22b) ^C(¯fi) =1hD∑d=1(⟨^Ci+12ed⟩−⟨^Ci−12ed⟩)=C(f)+O(ΔRq,Δvq∥,Δμq), (22c)

where are the orders of the schemes used to discretize the Vlasov and collision terms, respectively.

Equation (22) requires the computation of the face-averaged Vlasov and collisional fluxes defined by (20), for which we use the discretization described in collelaetal2011 . The following relationships between the cell-centered and face-centered quantities and cell-averaged and face-averaged quantities to fourth order will be used in the subsequent discussion:

 ⟨u⟩i±12ed =ui±12ed+h224D∑d′=1d′≠d∂2u∂ξ2d′∣∣ ∣∣i±12ed+O(h4), (23a) ¯ui =ui+h224D∑d=1∂2u∂ξ2d∣∣ ∣∣ξi+O(h4), (23b)

for an arbitrary variable , where denotes the face-averaged value, denotes the cell-averaged value, and and denote the face-centered and cell-centered values, respectively. We note that it is sufficient to compute the second derivatives in the RHS of (23) to second order using centered finite-differences since they are multiplied by . The face-averaged fluxes in (22) are computed to fourth order as

 ^Vdi±12ed= D∑s=1⟨Nsd⟩i±12ed⟨Vs⟩i±12ed +h2D∑s=1{G⊥,d0⟨Nsd⟩i±12ed}⋅{G⊥,d0⟨Vs⟩i±12ed}+O(h4), (24a) ^Cdi±12ed= D∑s=1⟨Nsd⟩i±12ed⟨Cs⟩i±12ed +h2D∑s=1{G⊥,d0⟨Nsd⟩i±12ed}⋅{G⊥,d0⟨Cs⟩i±12ed}+O(h4), (24b)

where and are the column vectors of the face-averaged metric quantities collelaetal2011 .

Equation (24a) requires computing . Writing each component of four-dimensional Vlasov term as an advection operator,

 (25)

where are the spatial dimensions corresponding to coordinates, respectively, the face-averaged flux is computed to fourth order from the discrete convolution as

 (26)

The face-averaged advective coefficients are computed by evaluating and at the face centers using (2) and transforming them to face-averaged quantities using (23a). The cell-averaged distribution function in the computational space is defined and related to as follows:

 ¯¯fi=1Vi∫\displaylimitsViB∗∥fdξ   ⇒¯¯fi=J−1i[¯fi−h212∇ξ¯f⋅∇ξJ]+O(h4). (27)

The derivatives in (26) and in (27) are computed using second-order central finite-differences since they are multiplied by . Finally, the face-averaged distribution function in the computational space is computed from cell-averaged values using the fifth-order weighted essentially nonoscillatory (WENO) scheme jiangshu with upwinding based on the sign of . Thus, the algorithm to compute in (22a) is fourth-order accurate, and in (22b).

Equation (24b) requires the computation of the face-averaged collision flux . We note that the velocity space grid is Cartesian in our drift-kinetic model and independent of the configuration space grid; thus, if , and (24b) only needs . The face-averaged flux is obtained from the face-centered collision flux

 Cdi±12ed=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩Γv∥i±12ev∥,d=dv∥Γμi±12eμ,d=dμ0,otherwise;   Γv∥i±12ev∥=→Γi±12ev∥⋅^v∥,Γμi±12eμ=→Γi±12eμ⋅^μ (28)

using (23a). Equation (4a) involves derivatives in the velocity space only, and to simplify the subsequent discussion, two-dimensional grid indices will be used, such that

 uk+p,l+q≡ui+pev∥+qeμ;  k=i⋅ev∥,  l=i⋅eμ (29)

for an arbitrary grid variable and . The computation of is described in the following paragraphs and can be extended trivially to computing .

Evaluating (4b) at the face , the Fokker-Planck collision flux along the dimension is given by

 Γv∥k+12,l≡Γv∥i+12ev∥=[σv∥f+κv∥v∥∂f∂v∥+κv∥μ∂f∂μ]k+12,l. (30)

The face-centered advective term is computed as

 [σv∥f]k+12,l=σv∥,k+12,lfk+12,l, (31)

where the distribution function at the interface is computed using the fifth-order upwind interpolation method based on the sign of , i.e.,

 σv∥,k+12,l <0, fk+12,l =130fk−2,l−1360fk−1,l+4760fk,l+2760fk+1,l−120fk+2,l, (32a) σv∥,k+12,l ≥0, fk+12,l =130fk+3,l−1360fk+2,l+4760fk+1,l+2760fk,l−120fk−1,l. (32b)

Note that the advective term is on the RHS of (1), and a negative value of implies a right-moving wave. The cell-centered values of the distribution function are obtained from its cell-averaged values using (23b). The diffusion terms are computed as follows:

 [κv∥v∥∂f∂v∥]k+12,l=κv∥v∥,k+12,l∂f∂v∥∣∣∣k+12,l, (33)

and

 [κv∥μ∂f∂μ]k+12,l=12⎛⎝κv∥μ,k+12,l−12∂f∂μ∣∣∣k+12,l−12+κv∥μ,k+12,l+12∂f∂μ∣∣∣k+12,l+12⎞⎠, (34)

where

 ∂f∂μ∣∣∣k+12,l±12=12⎛⎝∂f∂μ∣∣∣k,l±12+∂f∂μ∣∣∣k+1,l±12⎞⎠. (35)

The derivatives at the cell faces in (33) and (35) are computed using fourth-order central differences:

 ∂f∂v∥∣∣∣k+12,l=1Δv∥(124fk−1,l−98fk,l+98fk+1,l−124fk+2,l), (36)

and similarly for . The advective and diffusive coefficients in (31), (33), and (34) are related to the Rosenbluth potentials through (6). Appendix A describes the computation of the Rosenbluth potentials from the distribution function by discretizing and solving (7), which is defined on the infinite velocity domain but is solved on a truncated numerical domain. Our current implementation is second-order accurate, and therefore, the coefficients in (31), (33), and (34) are calculated by discretizing (6) with second-order finite differences:

 σv∥,k+12,l= 1Δv∥(−φk,l+φk+1,l), (37a) κv∥v∥,k+12,l= −12⎡⎣ϱk−1,l−2ϱk,l+ϱk+1,lΔv2∥+ϱk,l−2ϱk+1,l+ϱk+2,lΔv2∥⎤⎦, (37b) κv∥μ,k+12,l+12= −⎛⎝2μl+12mB⎞⎠ϱk,l−ϱk+1,l−ϱk,l+1+ϱk+1,l+1Δv∥Δμ. (37c)

The magnetic field magnitude in (37c) varies in the configuration space only. Overall, the algorithm to compute in (22a) is second-order accurate and in (22c).

#### Mass and Energy Conservation

The Fokker-Planck collision term conserves mass, momentum, and energy in its analytical form, and it is important that these quantities are discretely conserved as well. Mass conservation is expressed as

 ∫\displaylimitsΩvc(f)dv=0⇒∫\displaylimits∂Ωv→Γ⋅ds=0 (38)

where is the velocity space domain, and is the velocity space coordinate. Equation (38) is satisfied if the normal flux is zero over the entire velocity space boundary, and this is enforced by setting the coefficients , , , , and at the velocity boundaries to zero. Energy conservation is expressed as

 E{c(f)}(R)≡∫\displaylimitsΩvc(f)(v2∥+μB)dv=0, (39)

and is enforced using the approach described in taitanoetal2015 , adapted for Cartesian velocity coordinates. Let

 c(a)(f)=νc∇(v∥,μ)⋅[→σf],   c(d)(f)=νc∇(v∥,μ)⋅[↔κ∇(v∥,μ)f], (40)

and let and be the corresponding spatially discretized terms. Equation (39) implies

 E{c(d)(f)}(R)=−E{c(a)(f)}(R). (41)

We define

 νE(R)=−E{^C(a)(¯f)}(R)/E{^C(d)(¯f)}(R), (42)

and modify (4) as

 c(f)=νc∇(v∥,μ)⋅→Γ,     →Γ =→σf+νE(R)↔κ∇(v∥,μ)f (43)

We note that as for a consistent discretization of the collision term, and this factor counteracts the energy imbalance in the collision term resulting from the discretization errors. Equation (43) is discretized as decribed in the preceding discussions. Discrete enforcement of momentum conservation is also discussed in taitanoetal2015 ; however, the numerical test cases considered in this paper have negligible average velocities and yielded stable and convergent results without momentum conservation enforcement. In our implementation, discrete energy conservation is necessary to preserve a Maxwellian solution; without the modification described by (43), unphysical cooling of the Maxwellian solution is observed that is alleviated by refining the grid in the velocity space. This behavior is similar to that reported in taitanoetal2015 .

## 4 Time Integration

Equation (22a) represents a system of ODEs in time that are solved using high-order, multistage explicit and IMEX time integrators in COGENT. The plasma density and temperature vary significantly in the edge region along the radial direction, and as a result, the stiffness of (22a

) changes substantially from the inner edge region (adjacent to the core) to the outer edge. We use a simple problem to analyze the eigenvalues of the RHS of (

22a). The two-dimensional (in space) Cartesian domain is with specified density and temperature and a Maxwellian distribution function. The periodic domain is , where are the non-dimensional spatial coordinates; are the physical coordinates; and is the reference length. A uniform density is specified, and the temperature is specified as

 T=1+0.1cos(2πy), (44)

where and are nondimensional quantities. The magnetic field is