Low frequency fully kinetic simulation of the toroidal ion temperature gradient instability

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.¹ 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.² 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.² 

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%.⁴ In some H-mode edge pedestal, $ρi/L∼0.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,⁷ $ΩiΔt=0.2$⁠, and for microtearing simulations,⁸ $Ω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.⁹ FK ions/fluid electrons have been proposed for an extended-MHD closure,¹⁰ and have been used extensively to study kinetic-MHD stability in reversed field configurations.¹¹ A FK ion, gyrokinetic electron simulation model was first developed to study magnetic reconnection in space plasmas,¹² and FK simulations are used for studying RF waves in tokamaks.¹³ A first-order implicit δf FK ion drift-kinetic electron model was developed,² 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 model¹⁴ was later developed. This model was used to study the slab ITG instability.¹⁵ 

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.²² 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 f_i, given as

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

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

over all $x$ and $v$⁠, and the majority of the mass of f₀ 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 R_x, represented by a series expansion in ϵ as:

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

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

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

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

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 f₀, the perturbed part of the distribution function evolves as

A local Maxwellian is selected for f₀, where temperature and density are allowed to vary with R_x, that is

where T_i 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

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 T_i and $Ni$⁠, assuming constant values for the terms depending on R_x in Eqs. (8) and (9). We use the notation

so that $κT,κN>0$ correspond to equilibrium quantities that decrease with R_x. In addition, we now take T_i 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

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

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

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

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

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 = R₀. 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

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

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 r₀ as

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 r₀, the field-line-following coordinates are defined as:

where

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:

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

Fixed boundary conditions are taken in x as

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

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

where

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

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

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 N_p computational particles have been loaded according to a distribution function, f_L, selected to be an exact solution of Eq. (2). In this work, we have taken a spatially uniform Maxwellian distribution for f_L. The computational particles are taken to evolve according to the characteristics of Eq. (10) as

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

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

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

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

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 model^22,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:

where $R$ is the guiding center position, μ is the magnetic moment defined by $μ=v⊥2/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

In terms of cylindrical coordinates, this gives

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

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.

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

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

from which we arrive at

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ζ$⁠.

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×10−3, 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.

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 $ŝ≡r0q0q′0, LT=κT−1$⁠, and $LN=κN−1$⁠. In Fig. 5, a scan is performed measuring the real frequencies and growth rates as a function of the parameter $R0/LT$⁠, where $LT=κT−1$⁠. 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×10−3$ and $γ/Ωi=6.0×10−4$⁠, and for the global GEM gyrokinetic code, we have $ωR/Ωi=−2.26×10−3$ and $γ/Ωi=7.0×10−4$⁠. For the simulations in Fig. 5, the following resolution was used in the FK code: N_x = 128, N_y = 32, N_z = 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 N_x = 256, N_y = 64, and N_z = 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.

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

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

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

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

For our purposes, however, it is more convenient to have an integration scheme to be in terms of the cylindrical velocity components v_R, v_Z, and $vζ$⁠. A natural definition of the cylindrical velocity components in terms of p_R, p_Z, and $pζ$ comes from the continuous Lagrangian in Eq. (A1) by

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

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

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

and

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

yielding

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 $ζ̇¯$_:

Step 2: Advance coordinates via

Step 3: Advance velocity components via

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

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

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