A fully kinetic ion model is useful for the verification of gyrokinetic turbulence simulations in certain regimes, where the gyrokinetic model may break down due to the lack of small ordering parameters. However, for a fully kinetic model to be of value, it must first be able to accurately simulate low frequency drift-type instabilities typically well within the domain of gyrokinetics. Here, a fully kinetic ion model is formulated with weak gradient drive terms and applied to the toroidal ion-temperature-gradient (ITG) instability for the first time. Implementation in toroidal geometry is discussed, where orthogonal coordinates are used for particle dynamics, but field-line-following coordinates are used for the field equation allowing for high resolution of the field-aligned mode structure. Variational methods are formulated for integrating the equation of motion allowing for accuracy at a modest time-step size. Linear results are reported for both the slab and toroidal ITG instabilities. Good agreement with full Vlasov and gyrokinetic theory is demonstrated in slab geometry. Good agreement with global gyrokinetic simulation is also shown in toroidal geometry.

Gyrokinetic (GK) simulation models have been highly successful for modeling core microturbulence. A fully kinetic (FK) ion model that retains the cyclotron motion is promising because it would not be subject to the small expansion parameters used in the gyrokinetic ordering. On the other hand, if one is studying low frequency phenomena relative to the ion cyclotron frequency, the plethora of high frequency waves can make fully kinetic models intractable. Recently, a multiscale implicit method has been developed that captures the low frequency physics of the ion acoustic wave.1 A key challenge of this work is to avoid noise and stability constraints imposed by the ion Bernstein waves. Additionally, in an electromagnetic model, the compressional Alfvén wave will be present and will require an implicit method.2 In this paper, we extend the model in Ref. 1 to study the toroidal ion temperature gradient (ITG) instability. This is a problem where gyrokinetics performs extremely well and provides a significant challenge for a fully kinetic model due to the presence of high frequency ion Bernstein waves. The toroidal ITG instability provides a test bed to show the viability of a fully kinetic model. Such a simulation model can be extended to include electromagnetics and drift or gyrokinetic electrons.2 

A fully kinetic model is useful for the verification of gyrokinetics in situations where order parameters are not so small. There has been recent work on developing simulations using drift-cyclotron kinetic equations to study the breakdown of gyrokinetics using a slab four-dimensional continuum model with one of the variables being the gyrophase.3,4 It was found that the gyrokinetic approximation breaks down when density perturbations exceed 20%.4 In some H-mode edge pedestal, ρi/L0.2, where ρi is the ion gyroradius and L is the characteristic profile gradient scale.5,6 Additionally, for many topical problems in gyrokinetic simulation research, the required timestep turns out to be quite small. For example, for well-resolved gyrokinetic simulation of the DIII-D H-mode pedestal, ΩiΔt=1,5,6 where Ωi is the ion cyclotron frequency. For the simulation of NSTX core plasmas,7ΩiΔt=0.2, and for microtearing simulations,8ΩiΔt=0.25. A value of ΩiΔt=0.25 corresponds to the cyclotron period being approximately 25 times the timestep.

There is a long history of models using fully kinetic ions, also called “Lorentz ions,” for various applications.9 FK ions/fluid electrons have been proposed for an extended-MHD closure,10 and have been used extensively to study kinetic-MHD stability in reversed field configurations.11 A FK ion, gyrokinetic electron simulation model was first developed to study magnetic reconnection in space plasmas,12 and FK simulations are used for studying RF waves in tokamaks.13 A first-order implicit δf FK ion drift-kinetic electron model was developed,2 specifically for the purpose of studying low frequency microturbulence in magnetic fusion plasmas where a first-order implicit method was used to eliminate the constraint on the time step imposed by the fast compressional wave. A more accurate second-order implicit FK ion hybrid model14 was later developed. This model was used to study the slab ITG instability.15 

The toroidal ITG instability is a very well-studied problem in the context of gyrokinetics,16–21 and therefore provides a good starting point for the development of FK ion models which can be used to verify gyrokinetics. We develop a kinetic equilibrium model to include temperature and density gradients. To ensure the long time accuracy of the particle orbits in the FK ion model with a reasonable timestep size, we develop an integration scheme based on variational principles. Here, we focus on the “Cyclone DIII-D base case parameter set”17,21 as a benchmark for the FK toroidal ITG simulation model. Comparisons are made with both theory and the global GEM code.22 This paper is organized as follows, we begin by presenting the FK equations for modeling the ITG instability in Sec. II. In Sec. III, we discuss how the tokamak geometry is modeled using FK ions. In Sec. IV, the simulation model is discussed, and in Sec. V, the particle trajectories are compared directly with guiding center equations. Finally, Sec. VI presents the linear results in both slab and toroidal geometries as well as comparisons with theory and global gyrokinetic simulation.

In this section, we present the governing equations for the FK ion model of toroidal ITG instability. First, an approximate equilibrium distribution function, including temperature and density gradients, is constructed for the FK ion model, assuming equilibrium quantities to vary on length scales which are much larger than a typical ion gyro-radius. Our starting point is the Vlasov equation for the evolution of the ion distribution function fi, given as

fit+v·fi+qimi(E+v×B)·vfi=0,
(1)

where qi and mi are the ion charge and mass, respectively, E is the electric field, and B is the magnetic field. It is assumed that fi can be separated into equilibrium and perturbed parts as fi=f0+δf. For consistency with Eq. (1), f0 is required to satisfy

v·f0+qimi(v×B)·vf0=0.
(2)

We seek an approximate equilibrium distribution function f0ϵ, such that f0ϵf0 for some f0 satisfying Eq. (2), when ϵ1. To define the small parameter ϵ, we assume reference quantities Leq,Bref, and vref, such that

Leqf0f01,BBref1,

over all x and v, and the majority of the mass of f0 is contained in |v|vref. In terms of these quantities, we define ϵ=ρref/Leq, with ρref=mivref/qiBref. We note that the first term on the left hand side of Eq. (2) is O(ϵ) smaller than the second term, following directly from the definition of ϵ. Next, we take a curvilinear coordinate, x, such that x·B=0. In our numerical implementation, we will choose x to be the radial filed-line-following coordinate introduced in Sec. III. The approximate equilibrium distribution function is then taken to depend on the quantity Rx, represented by a series expansion in ϵ as:

Rx=n=0ϵnr(n),
(3)

where r(0)=x, and the remaining coefficients satisfy

(v×B)·vr(n+1)=BrefLeqvrefv·r(n).

We take RxN to be the series in Eq. (3) truncated past O(ϵN). Our approximate equilibrium distribution function is then f0ϵ(RxN), which can be shown to satisfy

v·f0ϵ+qimi(v×B)·vf0ϵ=ϵNv·r(N)f0ϵRxN.

Heuristically, from this expression, we expect Eq. (2) to be satisfied by f0ϵ(RxN) in the limit ϵ0+ or in the limit N for small, but finite ϵ.

For the remainder of the paper, we take the expansion in Eq. (3) to first order in ϵ, which can be written as

Rx=x+miqiv×b̂B·x,
(4)

where B is the magnitude of B,b̂ is the unit vector field in the direction of B, and we have dropped the superscript on the left hand side of Eq. (4) for the simplicity of notation. This first order approximation is interpreted as the projection of the particle coordinate x to the guiding center coordinate, and an examination of the terms on the right hand side of Eq. (4) shows the second term to be O(ϵ) smaller than the first term. To construct the equilibrium distribution function, we also consider the ion kinetic energy, given by

K=mi2v2,
(5)

where v2=v·v. It is straightforward to show that any function depending only on K will satisfy Eq. (2) exactly. An approximate equilibrium distribution function is then constructed from the quantities in Eqs. (4) and (5) as f0=f0(Rx,K), which satisfies Eq. (2) to first order in ϵ. Here, we have dropped the superscript ϵ in the approximate equilibrium distribution function for the simplicity of notation. With this form for f0, the perturbed part of the distribution function evolves as

δft+v·δf+qimi(E+v×B)·vδf=qiE·vf0Kx·(E×b̂B)f0Rx.
(6)

A local Maxwellian is selected for f0, where temperature and density are allowed to vary with Rx, that is

f0(K,Rx)=Ni(Rx)(2πTi(Rx)/mi)3/2eK/Ti(Rx),
(7)

where Ti and Ni are the ion equilibrium temperature and the density, respectively. With this form, the partial derivatives on the right hand side of Eq. (6) are given by

f0K=1Tif0,
(8)
f0Rx=(NiRx/Ni+(KTi32)TiRx/Ti)f0.
(9)

In the present work, we consider only the linearized version of Eq. (6), that is, we neglect the term containing the electric field on the left hand side to study stability in the presence of a small perturbation. In addition, we do not consider the effects due to profile variations of Ti and Ni, assuming constant values for the terms depending on Rx in Eqs. (8) and (9). We use the notation

κT=TiRx/Ti|Rx=0,κN=NiRx/Ni|Rx=0,

so that κT,κN>0 correspond to equilibrium quantities that decrease with Rx. In addition, we now take Ti and Ni to represent the constant values Ti(Rx=0) and Ni(Rx=0), respectively. With this notation, the linearized version of Eq. (6), assuming Eq. (7) and neglecting equilibrium profile variations, is written as

δft+v·δf+qimi(v×B)·vδf=qiTiE·vf0+(E×b̂)·xB[κN+(miv22Ti32)κT]f0.
(10)

To close Eq. (10), an electrostatic field is assumed as

E=ϕ,
(11)

where ϕ is the electrostatic potential. Assuming adiabatic electrons and quasineutrality, an equation for ϕ is given as

eϕTe=δnNi,
(12)

where e and Te are the electron charge and the temperature, respectively, and δn is obtained by

δn=3δfd3v.

In this section, we first introduce the basic toroidal coordinate system used to express the magnetic field model for this study. A field-line-following coordinate system is then defined in terms of basic toroidal coordinates, which are used to express the so-called flux tube domain over which our model equations are solved.16,17,22 Since the toroidal ITG mode is characterized by long wavelengths parallel to the magnetic field and short perpendicular wavelengths, a coordinate system which is aligned with the magnetic field allows for coarser resolution in the direction of the magnetic field. In addition, it allows us to define a minimal simulation volume necessary to capture the relevant physics.

The basic toroidal coordinate system is most easily defined in terms of a cylindrical coordinate system (R,Z,ζ).23,24 Here, R is the cylindrical radius and Z is the cylindrical axis. The angle ζ is taken to be the negative azimuthal angle, i.e., increasing ζ corresponds to a clockwise rotation about the Z axis, so that (R,Z,ζ) is right handed. The Cartesian coordinates (x1,x2,x3) are expressed in terms of (R,Z,ζ) by

x1=Rcosζ,
x2=Rsinζ,
x3=Z.

In the basic toroidal coordinate system, we refer to the Z axis of the cylindrical system as the major axis and ζ as the toroidal angle. The minor axis is a circle in the Z = 0 plane with radius R = R0. The toroidal coordinates consist of the toroidal angle ζ, the minor radius r measuring the distance from the minor axis in a constant ζ plane, and an angle about the minor axis θ called the poloidal angle. The cylindrical coordinates R and Z are expressed in terms of r and θ by

R=R0+rcosθ,
Z=rsinθ.

Written in mixed cylindrical and toroidal coordinates, the assumed form of the magnetic field in our model is given by

B=B0(R0Rζ̂+rq(r)Rθ̂),
(13)

where the function q(r) is the safety factor. We take the form for q(r) to be a local expansion about a reference minor radius r0 as

q(r)=q0+(rr0)q0(r0).

The simulation domain is defined in terms of a field-line-following coordinate system. In terms of the toroidal coordinates and a reference minor radius r0, the field-line-following coordinates are defined as:

x=rr0,
(14)
y=r0q0(0θq̂(r,θ)dθζ),
(15)
z=q0R0θ,
(16)

where

q̂(r,θ)=ζ·Bθ·B.

These coordinates are non-orthogonal and are chosen to have the properties B·x=B·y=0. The simulation domain is taken to be a rectangular region in the field-line-following coordinates given as:

D={(x,y,z)|lx2xlx2,ly2yly2,q0R0πzq0R0π},

which corresponds to a long, thin tube in physical space following the twisted magnetic field.

The boundary conditions for a perturbed quantity A over the flux tube domain are as follows. Periodicity is enforced in y for fixed x and z as

A(x,y+ly,z)=A(x,y,z).
(17)

Fixed boundary conditions are taken in x as

A(lx2,y,z)=A(lx2,y,z)=0.
(18)

The boundary condition at the end points in z requires some discussion. In selecting the boundary condition, we wish to enforce the periodicity in θ for a fixed toroidal angle. In particular, we require

A(x(r),y(r,θ+2π,ζ),z(θ+2π))=A(x(r),y(r,θ,ζ),z(θ)).

From the coordinate transformation in Eqs. (14)–(16), we have

y(r,θ+2π,ζ)=y(r,θ,ζ)+δy(r),
z(θ+2π)=z(θ)+2πq0R0,

where

δy(r)=r0q0ππq̂(r,θ)dθ.

The boundary condition in z can therefore be written as the periodicity constraint

A(x,y+δy(x),z+lz)=A(x,y,z),
(19)

where lz=2πq0R0. An example of a function satisfying Eqs. (17)–(19) can be given by

A(x,y,z)=C(x)ei2πnly(yδy(x)zlz)ei2πmlzz,

where C(lx2)=C(lx2)=0.

In this section, we present the details of numerical implementation of our model. We employ a δf particle-in-cell (PIC) method to solve Eq. (10) for the perturbed ion distribution function.25–29 Although it is trivial to obtain ϕ from Eq. (12), some care must be taken to account for the toroidal geometry when computing the gradient in Eq. (11). In addition, the toroidal geometry requires modification in the particle loading step and in depositing the particles to the grid. Details of the integration scheme applied to the particle equations of motion are given in Sec. V and the  Appendix.

In our simulations, we employ the δf PIC method to evolve the perturbed part of the ion distribution function in Eq. (10). We assume that Np computational particles have been loaded according to a distribution function, fL, selected to be an exact solution of Eq. (2). In this work, we have taken a spatially uniform Maxwellian distribution for fL. The computational particles are taken to evolve according to the characteristics of Eq. (10) as

dxpdt=vp,
(20)
dvpdt=qimi(vp×B(xp)),
(21)

for p=1,,Np. By choosing fL to satisfy Eq. (2) with each particle following Eqs. (20) and (21), the computational particles will be distributed according to fL for all times throughout the simulation. Next, weights are defined for each computational particle as

wp=δf(xp,vp)fL(xp,vp).
(22)

By taking the time derivative of Eq. (22), we obtain an equation for the evolution of wp from Eqs. (20) and (21) and Eq. (10) as

dwpdt=(qiTiE·v+(E×b̂)·xB[κN+(miv22Ti32)κT])f0fL|xp,vp.
(23)

In our implementation, cylindrical coordinates are used to evolve the particles in phase space according to Eqs. (20) and (21). Particle positions are transformed back to the field-line-following coordinates to be deposited to the grid and to interpolate the field values in Eq. (23).

The computational particles evolve along the equilibrium trajectories according to Eqs. (20) and (21). The choice of coordinates used to describe a particle's position will determine how simply Eqs. (20) and (21) can be expressed. Since our computational domain uses a uniform grid in the field-line-following coordinates, the deposit stage requires particle positions to be in the field-line-following coordinates, as well. Expressing Eqs. (20) and (21) in terms of (x, y, z), however, results in a cumbersome form involving complicated geometric terms. In this work, we have made the choice to use the cylindrical (R,Z,ζ) coordinates in the particle push stage and use the analytical coordinate transformations of Eqs. (14)–(16) to convert to (x, y, z) in the deposit stage. In terms of cylindrical coordinates, Eqs. (20) and (21) can be written as

R¨=Rζ̇2+qimi(ŻBζRζ̇BZ),
(24)
Z¨=qimi(Rζ̇BRṘBζ),
(25)
Rζ¨=2Ṙζ̇+qimi(ṘBZŻBR),
(26)

where dots are used to denote a derivative with respect to time, and we have assumed B to be expressed as

B=BRR̂+BZẐ+Bζζ̂.

Using cylindrical coordinates in the particle push has an additional advantage if more complicated magnetic geometries are to be used. For example, if the more general Miller equilibrium model22,30 is used, the integration scheme would not need to be changed. The only changes needed would be in the transformation to (x, y, z) and the expression of B in the cylindrical coordinates. A discrete variational integration scheme is applied to the fully kinetic equations of motion Eqs. (24)–(26), and the details are presented in the  Appendix.

To test the fully kinetic, variational integration scheme, we consider the motion of a trapped particle in a tokamak field. For comparison, we also consider a predictor-corrector, guiding the center particle integrator to compute the trapped particle orbit following the equations of motion:

dRdt=vb̂+miqiBb̂×[v2(b̂·R)b̂+μBlnBR],dvdt=v22(b̂·lnBR),

where R is the guiding center position, μ is the magnetic moment defined by μ=v2/2B with v being the magnitude of velocity perpendicular to b̂ and v the velocity component along b̂. In Table I, the initial conditions for both the fully kinetic and guiding center integrators are given. For consistency, the initial conditions of the guiding center integrator are related to the initial conditions of the fully kinetic integrator through

R=x+miqiv×b̂B.

In terms of cylindrical coordinates, this gives

RgcRfk+miqiv×b̂B·R,
ZgcZfk+miqiv×b̂B·Z,
ζgcζfk+miqiv×b̂B·ζ,

where the subscripts “gc” and “fk” indicate the guiding center and fully kinetic coordinates, respectively. We use the magnetic field model in Eq. (13) expressed in cylindrical components as

B=B0R0Rq(r)(ZR0R̂+(RR0)R0Ẑ+q(r)ζ̂).
(27)

The parameters of the magnetic field model are the same as in Table II. For both integrators, a time step size of ΩiΔt=0.125 is taken, and the simulations are run to Ωit=2.5×104. The orbits in the (R, Z) plane along with the full 3D orbits mapped to the Cartesian coordinates (x1,x2,x3) are shown in Fig. 1. Both integrators correctly capture the trapped particle orbit and agree well with each other.

TABLE I.

Initial conditions for the trapped particle orbit.

CoordinateFully kineticGuiding center
R/R0 1.176 1.180 
Z/R0 3.912 ×103 0.000 
ζ −4.258 ×104 0.000 
vR/R0Ωi 3.371 ×103 … 
vZ/R0Ωi 3.371 ×103 … 
vζ/R0Ωi 1.798 ×103 … 
v/R0Ωi … 2.212 ×103 
B0μ/(R0Ωi)2 … 1.232 ×105 
CoordinateFully kineticGuiding center
R/R0 1.176 1.180 
Z/R0 3.912 ×103 0.000 
ζ −4.258 ×104 0.000 
vR/R0Ωi 3.371 ×103 … 
vZ/R0Ωi 3.371 ×103 … 
vζ/R0Ωi 1.798 ×103 … 
v/R0Ωi … 2.212 ×103 
B0μ/(R0Ωi)2 … 1.232 ×105 
TABLE II.

Cyclone DIII-D base case parameter set.

R0/ρir0/R0q0ŝR0/LTR0/LNτ
445.0 0.18 1.4 0.78 6.9 2.2 1.0 
R0/ρir0/R0q0ŝR0/LTR0/LNτ
445.0 0.18 1.4 0.78 6.9 2.2 1.0 
FIG. 1.

Comparison of a FK variational integrator with a guiding center integrator.

FIG. 1.

Comparison of a FK variational integrator with a guiding center integrator.

Close modal

The conservation properties of our fully kinetic integrator are tested with two exact constants of motion that can be found for the system Eqs. (20) and (21), using the form of B from Eq. (13). These are the kinetic energy K defined in Eq. (5) and the toroidal angular momentum pζ defined in Eq. (A8). The kinetic energy can be expressed in terms of cylindrical velocity components simply as

K=mi2(vR2+vZ2+vζ2).

The toroidal angular momentum is conserved due to the symmetry of B with respect to the toroidal angle. Its computation, however, requires the toroidal component of a vector potential. From Eqs. (A15)–(A17), it can be verified that a vector potential for the field model in Eq. (27) is given by

A=B0R0ZRR̂B0Rr0rrq(r)drζ̂,

from which we arrive at

pζ=miRvζqiB0r0rrq(r)dr.

In Fig. 2, we examine the relative variation of K and pζ over a time period of 2.5×104 for different time step sizes. The initial conditions are taken from Table I. The amount by which the conservation of these quantities is violated is characterized by bounded oscillations, with amplitudes that decrease with the time step size. Figure 3 shows the amplitude of these oscillations as a function of the time step size. The amplitudes are shown to decrease as O(Δt2) for both K and pζ.

FIG. 2.

Conservation of kinetic energy and toroidal angular momentum in a fully kinetic integrator.

FIG. 2.

Conservation of kinetic energy and toroidal angular momentum in a fully kinetic integrator.

Close modal
FIG. 3.

Convergence of kinetic energy and toroidal momentum conservation with ΩiΔt.

FIG. 3.

Convergence of kinetic energy and toroidal momentum conservation with ΩiΔt.

Close modal

As a first test, we simulate the linear ITG instability in slab geometry using the FK model equations. This is a useful benchmark since the derivation of linear FK and GK dispersion relations is tractable. These dispersion relations can be used to compare simulation results against. We use the test case from Ref. 15 with parameters τ=4.0,kρi=2.0×103,kρi=0.2,κN=0.0 and scan over the values of κTρi. Here, τ is defined as τ=qiTe/eTi. Figure 4 shows the results of simulation compared with the solutions to both FK and GK dispersion relations. The simulations were performed using Δy=0.65ρi,Δz=65.5ρi,ΩiΔt=0.125 and 4 194 304 particles. Excellent agreement is observed between the FK simulations, the FK dispersion solutions, and the GK dispersion relations.

FIG. 4.

Simulation results for the slab ITG FK ion model compared to the dispersion theories from the FK ion model and the GK ion model.

FIG. 4.

Simulation results for the slab ITG FK ion model compared to the dispersion theories from the FK ion model and the GK ion model.

Close modal

Next, we present the simulation results for the FK ion model applied to the ITG instability model. The “Cyclone DIII-D base case parameter set,” given in Table II, is used to benchmark our algorithm, where ŝr0q0q0,LT=κT1, and LN=κN1. In Fig. 5, a scan is performed measuring the real frequencies and growth rates as a function of the parameter R0/LT, where LT=κT1. The domain size is taken to be lx=32ρi,ly=12.57ρi,lz=3914.4ρi, and all other parameters are taken from Table II. Comparisons are made with the global GEM code and good agreement is observed. For the case with R0/LT=6.9, the real frequencies of the two codes agree within 6 percent and the growth rates agree within 17 percent. For the FK code, we have ωR/Ωi=2.40×103 and γ/Ωi=6.0×104, and for the global GEM gyrokinetic code, we have ωR/Ωi=2.26×103 and γ/Ωi=7.0×104. For the simulations in Fig. 5, the following resolution was used in the FK code: Nx = 128, Ny = 32, Nz = 48, corresponding to Δx=0.24ρi,Δy=0.39ρi, and Δz=81.55ρi. A time step size of ΩiΔt=0.2 was used with 21 particles per cell. An additional test was performed for the case of R0/LT=6.9 using the resolution Nx = 256, Ny = 64, and Nz = 96 with ΩiΔt=0.125 and 21 particles per cell. This produced no significant difference in the real frequency or damping rate, suggesting that our FK ion model implementation is converged for the results in Fig. 5.

FIG. 5.

Comparisons of fully kinetic simulation with global gyrokinetic (GEM) simulations. Real frequencies and growth rates are shown as a function of R0/LT.

FIG. 5.

Comparisons of fully kinetic simulation with global gyrokinetic (GEM) simulations. Real frequencies and growth rates are shown as a function of R0/LT.

Close modal

In this paper, a fully kinetic ion model was developed to simulate low-frequency microinstabilities in toroidal geometry. Key to this work was the inclusion of weak equilibrium gradients in the FK ion model as well as an inhomogeneous toroidal magnetic field, as was presented in Sec. II. The full geometric effects in toroidal flux tube geometry have been accounted for in the FK implementation. To make a FK simulation viable, a particle integration scheme has been developed, based on variational principles, which is shown to provide accurate and stable orbits over long simulation times. Additionally, the integration scheme has been shown to accurately conserve the constants of motion along a particle's trajectory. FK ion simulations of toroidal ITG instability have been carried out for the first time, and good agreement obtained between the FK ion implementation and the global GEM gyrokinetic code. Slab comparisons with theory show good agreement as well. This work provides a starting point for developing a FK ion model that can be used to verify gyrokinetic codes. Future work includes exploring the practical limits of gyrokinetics, as well as adding more realism to the FK simulation model, such as including nonlinearities, electromagnetic perturbations, and kinetic electrons.

Research was supported by the U.S. Department of Energy.

Rather than a direct discretization of Eqs. (24)–(26), e.g., via Runge-Kutta, we derive a scheme based on discrete variational methods.31–34 Variational integrators have been developed for guiding center equations in magnetic fusion applications.35,36 Integration schemes derived from variational methods automatically have a number of desirable properties. In particular, when a strict variational formalism is followed, the resulting algorithm will be symplectic; it will exactly preserve the momentum associated with symmetries of the system and will exhibit excellent energy stability over long simulation times. Integration schemes based on variational principles are a broad area of research, which we do not attempt to summarize here. The starting point for our integrator is the position-momentum form of discrete Euler-Lagrange equations, for which a derivation can be found in Ref. 31. The Lagrangian for the single particle motion described in Eqs. (24)–(26) is written in the cylindrical coordinates as

L=mi2(Ṙ2+Ż2+R2ζ̇2)+qi(ṘAR+ŻAZ+Rζ̇Aζ),
(A1)

where A is the vector potential. It is noted that an integration scheme for the full nonlinear form of Eq. (6) could be derived along the same lines presented here by including the contribution of qiϕ to the Lagrangian. To form the discrete Lagrangian, we use the trapezoidal rule along with the approximations

ṘṘ¯=Rν+1RνΔt,
ŻŻ¯=Zν+1ZνΔt,
ζ̇ζ̇¯=ζν+1ζνΔt,

for t[νΔt,(ν+1)Δt]. This choice of quadrature will result in a scheme similar to the velocity Verlet method.37 The discrete Lagrangian over this interval is then

Ld=Δtmi2[Ṙ¯2+Ż¯2+(Rν)2+(Rν+1)22ζ̇¯2]+Δtqi2[Ṙ¯(ARν+ARν+1)+Ż¯(AZν+AZν+1)+ζ̇¯(RνAζν+Rν+1Aζν+1)],
(A2)

where Aν is understood to be the evaluation of A at the particle's location at time t=νΔt. From the position-momentum form of discrete Euler-Lagrange equations,31,33,34 an integration scheme in terms of discrete conjugate momenta can be given by

pRν=LdRν,pRν+1=LdRν+1,
(A3)
pZν=LdZν,pZν+1=LdZν+1,
(A4)
pζν=Ldζν,pζν+1=Ldζν+1.
(A5)

For our purposes, however, it is more convenient to have an integration scheme to be in terms of the cylindrical velocity components vR, vZ, and vζ. A natural definition of the cylindrical velocity components in terms of pR, pZ, and pζ comes from the continuous Lagrangian in Eq. (A1) by

pR=LṘ=miṘ+qiAR=mivR+qiAR,
(A6)
pZ=LŻ=miŻ+qiAZ=mivZ+qiAZ,
(A7)
pζ=Lζ̇=miR2ζ̇+qiRAζ=miRvζ+qiRAζ.
(A8)

The integration scheme in Eqs. (A3)–(A5) written in terms of velocity is then

vRν=qimiARν1miLdRν,
(A9)
vRν+1=qimiARν+1+1miLdRν+1,
(A10)
vZν=qimiAZν1miLdZν,
(A11)
vZν+1=qimiAZν+1+1miLdZν+1,
(A12)
vζν=qimiAζν1miRνLdζν,
(A13)
vζν+1=qimiAζν+1+1miRν+1Ldζν+1.
(A14)

Furthermore, it is convenient to be able to work directly with B, rather than having to first form a vector potential A. By expressing the curl of A in cylindrical coordinates, we have

BR=AζZ1RAZζ,
(A15)
BZ=1R(ARζR(RAζ)),
(A16)
Bζ=AZRARZ.
(A17)

Next, we illustrate how the components of A can be eliminated in Eqs. (A9) and (A10) in favor of the components of B by means of an approximation. Computing the partial derivatives of the discrete Lagrangian in Eq. (A2) gives

vRν=Ṙ¯Δt2Rν(ζ̇¯)2+qi2mi(ARν+1ARν)qiΔt2mi[Ṙ¯ARνRν+Ż¯AZνRν+ζ̇¯Rν(RνAζν)],
(A18)

and

vRν+1=Ṙ¯+Δt2Rν+1(ζ̇¯)2+qi2mi(ARνARν+1)+qiΔt2mi[Ṙ¯ARν+1Rν+1+Ż¯AZν+1Rν+1+ζ̇¯Rν+1(Rν+1Aζν+1)].
(A19)

In Eqs. (A18) and (A19), we make the following approximations, respectively:

ARν+1ARν+Δt(Ṙ¯ARνRν+Ż¯ARνZν+ζ̇¯ARνζν),ARνARν+1Δt(Ṙ¯ARν+1Rν+1+Ż¯ARν+1Zν+1+ζ̇¯ARν+1ζν+1),

yielding

vRν=Ṙ¯Δt2Rν(ζ̇¯)2qiΔt2mi[Ż¯Bζνζ̇¯RνBZν],vRν+1=Ṙ¯+Δt2Rν+1(ζ̇¯)2+qiΔt2mi[Ż¯Bζν+1ζ̇¯Rν+1BZν+1].

It is noted that this approximation does not follow the strict variational formalism, and hence we cannot guarantee that theories derived for variational integrators should hold when it is employed. Numerical results, from single particle motion, however, seem to exhibit excellent conservation properties for the constant of motion despite the approximation. Equations (A11)–(A14) can be approximated in the same manner to eliminate components of A, yielding the following procedure for computing (Rν+1,Zν+1,ζν+1) and (vRν+1,vZν+1,vζν+1) giving (Rν,Zν,ζν) and (vRν,vZν,vζν):

Step 1: Solve the following system of equations for the unknowns Ṙ¯,Ż¯, and ζ̇¯:

vRν=Ṙ¯ΔtRν2(ζ̇¯)2Δtqi2mi[Ż¯Bζνζ̇¯RνBZν],
(A20)
vZν=Ż¯+Δtqi2mi[Ṙ¯Bζνζ̇¯RνBRν],
(A21)
vζν=ζ̇¯Rν+Δtζ̇¯[Ṙ¯+Δt(Ṙ¯)22Rν]Δtqi2mi[Ṙ¯BZνŻ¯BRν].
(A22)

Step 2: Advance coordinates via

Rν+1=Rν+ΔtṘ¯,
(A23)
Zν+1=Zν+ΔtŻ¯,
(A24)
ζν+1=ζν+Δtζ̇¯.
(A25)

Step 3: Advance velocity components via

vRν+1=Ṙ¯+ΔtRν+12(ζ̇¯)2+Δtqi2mi[Ż¯Bζν+1ζ̇¯Rν+1BZν+1],
(A26)
vZν+1=Ż¯Δtqi2mi[Ṙ¯Bζν+1ζ̇¯Rν+1BRν+1],
(A27)
vζν+1=ζ̇¯(Rν)2+(Rν+1)22Rν+1+Δtqi2mi[Ṙ¯BZν+1Ż¯BRν+1].
(A28)

Note that this procedure requires the solution of a nonlinear system of three equations for the unknowns Ṙ¯,Ż¯, and ζ̇¯ in step 1. Here, this is accomplished by Picard iterations. Denoting the unknowns by x¯[Ṙ¯,Ż¯,ζ̇¯]T, the known velocity components by v¯[vRν,vZν,vζν]T, and taking the operator A(x¯) to represent the right hand side of Eqs. (A20)–(A22), we can express these equations as A(x¯)=v¯. We note that A can be decomposed into linear and nonlinear parts as

A(x¯)=Lx¯+(x¯).

In particular, consists of the second terms of Eqs. (A20) and (A22). The Picard iteration method proceeds as follows:

x¯0=L1v¯,x¯i=L1(v¯(x¯i1)),i=1,2,

We have found this procedure to converge rapidly, with two iterations being sufficient for a wide range of time step sizes.

1.
B. J.
Sturdevant
,
S. E.
Parker
,
Y.
Chen
, and
B. B.
Hause
,
J. Comput. Phys.
316
,
519
(
2016
).
2.
Y.
Chen
and
S. E.
Parker
,
Phys. Plasmas
16
,
052305
(
2009
).
3.
R. E.
Waltz
and
Z.
Deng
,
Phys. Plasmas
20
,
012507
(
2013
).
4.
Z.
Deng
and
R. E.
Waltz
,
Phys. Plasmas
22
,
056101
(
2015
).
5.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
,
Z.
Yan
,
R. J.
Groebner
, and
R. B.
Snyder
,
Phys. Rev. Lett.
109
,
185004
(
2012
).
6.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
,
R. J.
Groebner
,
Z.
Yan
,
A. Y.
Pankin
, and
S. E.
Kruger
,
Phys. Plasmas
20
,
055902
(
2013
).
7.
Y.
Chen
,
S. E.
Parker
,
G.
Rewoldt
,
S. H.
Ku
,
G. Y.
Park
, and
C. S.
Chang
,
Phys. Plasmas
15
,
055905
(
2008
).
8.
J.
Chowdhury
,
Y.
Chen
,
W.
Wan
,
S. E.
Parker
,
W.
Guttenfelder
, and
J. M.
Canik
,
Phys. Plasmas
23
,
012513
(
2016
).
9.
M.
Tanaka
,
J. Comput. Phys.
107
,
124
(
1993
).
10.
D.
Barnes
,
J.
Cheng
, and
S.
Parker
,
Phys. Plasmas
15
,
055702
(
2008
).
11.
E. V.
Belova
,
S. C.
Jardin
,
H.
Ji
,
M.
Yamada
, and
R.
Kulsrud
,
Phys. Plasmas
7
,
4996
(
2000
).
12.
Y.
Lin
,
Z.
Wang
,
Z.
Lin
, and
L.
Chen
,
Plasma Phys. Controlled Fusion
47
,
657
(
2005
).
13.
A.
Kuley
,
Z.
Lin
,
J.
Bao
,
X. S.
Wei
,
Y.
Xiao
,
W.
Zhang
,
G.
Sun
, and
N. J.
Fisch
,
Phys. Plasmas
22
,
102515
(
2015
).
14.
J.
Cheng
,
S. E.
Parker
,
Y.
Chen
, and
D.
Uzdensky
,
J. Comput. Phys.
245
,
364
(
2013
).
15.
D. D.
Schnack
,
J.
Cheng
,
D. C.
Barnes
, and
S. E.
Parker
,
Phys. Plasmas
20
,
062106
(
2013
).
16.
M. A.
Beer
, “
Gyrofluid models of turbulent transport in tokamaks
,” Ph.D. thesis,
Princeton University
,
1995
.
17.
S. E.
Parker
,
C.
Kim
, and
Y.
Chen
,
Phys. Plasmas
6
(
5
),
1709
(
1999
).
18.
S. E.
Parker
,
W. W.
Lee
, and
R.
Santoro
,
Phys. Rev. Lett.
71
,
2042
(
1993
).
19.
R.
Sydora
,
V.
Decyk
, and
J.
Dawson
,
Plasma Phys. Controlled Fusion
38
,
A281
(
1996
).
20.
S. E.
Parker
,
H.
Mynick
,
M.
Artun
,
V.
Decyk
,
J.
Kepner
,
W. W.
Lee
, and
W.
Tang
,
Phys. Plasmas
3
,
1959
(
1996
).
21.
A. M.
Dimits
,
G.
Bateman
,
M. A.
Beer
,
B. I.
Cohen
,
W.
Dorland
,
G. W.
Hammett
,
C.
Kim
,
J. E.
Kinsey
,
M.
Kotschenreuther
,
L. L.
Kritz
,
A. H.
Lao
,
J.
Mandrekas
,
W. M.
Nevins
,
S. E.
Parker
,
A. J.
Redd
,
D. E.
Shumaker
,
R.
Sydora
, and
J.
Weiland
,
Phys. Plasmas
7
(
3
),
969
(
2000
).
22.
Y.
Chen
and
S. E.
Parker
,
J. Comput. Phys.
220
,
839
855
(
2007
).
23.
R. D.
Hazeltine
and
J. D.
Meiss
,
Plasma Confinement
(
Dover Publications
,
2003
).
24.
W.
D'haeseleer
,
W.
Hitchon
,
J.
Callen
, and
J.
Shohet
,
Flux Coordinates and Magnetic Field Structure
(
Springer-Verlag
,
Berlin, Heidelberg
,
1991
).
25.
S. E.
Parker
and
W. W.
Lee
,
Phys. Fluids B
5
,
77
(
1993
).
26.
A. M.
Dimits
and
W. W.
Lee
,
J. Comput. Phys.
107
,
309
(
1993
).
27.
G.
Hu
and
J. A.
Krommes
,
Phys. Plasmas
1
,
863
(
1994
).
28.
A. Y.
Aydemir
,
Phys. Plasmas
1
,
822
(
1994
).
29.
M.
Kotschenreuther
,
Bull. Am. Phys. Soc.
34
,
2107
(
1993
).
30.
R. L.
Miller
,
M. S.
Chu
,
J. M.
Greene
,
Y. R.
Lin-Lui
, and
R. E.
Waltz
,
Phys. Plasmas
5
,
973
(
1998
).
31.
J. E.
Marsden
and
M.
West
,
Acta Numerica
10
,
357
(
2001
).
32.
E.
Hairer
,
C.
Lubich
, and
G.
Wanner
,
Geometric Numerical Integration
(
Springer-Verlag
,
Berlin, Heidelberg
,
2002
).
33.
M.
West
, “
Variational integrators
,” Ph.D. thesis,
California Institute of Technology
,
2004
.
34.
M.
Krauss
, “
Variational integrators in plasma physics
,” Ph.D. thesis,
Technishe Universitat Munchen
,
2013
.
35.
H.
Qin
and
X.
Guan
,
Phys. Rev. Lett.
100
,
035006
(
2008
).
36.
H.
Qin
,
X.
Guan
, and
W. M.
Tang
,
Phys. Plasmas
16
,
042510
(
2009
).
37.
L.
Verlet
,
Phys. Rev.
159
,
98
(
1967
).