Fully Kinetic Simulation of Ion-Temperature-Gradient Instabilities in Tokamaks

The feasibility of using full ion kinetics, instead of gyrokinetics, in simulating low-frequency Ion-Temperature-Gradient (ITG) instabilities in tokamaks has recently been demonstrated. The present work extends the full ion kinetics to the nonlinear regime and investigates the nonlinear saturation of a single-n ITG instability due to the E× B trapping mechanism (n is the toroidal mode number). The saturation amplitude predicted by the E× B trapping theory is found to agree with the saturation level observed in the simulation. In extending to the nonlinear regime, we developed a toroidal Boris full orbit integrator, which proved to be accurate in capturing both the short-time scale cyclotron motion and long time scale drift motion, with good kinetic energy conservation and toroidal angular momentum conservation in tokamak equilibrium magnetic fields. This work also extends the previous work from analytic circular magnetic equilibria to general numerical magnetic equilibria, enabling simulation of realistic equilibria reconstructed from tokamak experiments.


Introduction
Ion temperature gradients in tokamaks provide free energy to micro-instabilities called Ion-Temperature-Gradient (ITG) instabilities [1].The nonlinear development of these instabilities, i.e., ITG turbulence, is believed to play an important role in regulating particle and heat transport in tokamaks [2][3][4].There are numerous papers devoted to the gyrokinetic simulation of ITG turbulence, which employ gyrokinetic theory to decouple the high-frequency gyro-motion of ions from the low-frequency ITG modes [5][6][7][8][9][10].However, gyrokinetics rely on ordering assumptions in deriving the gyrokinetic equation.One of these ordering assumptions, ρ i /L n 1, becomes questionable in tokamak edge with steep density profile, where ρ i is the gyro-radius of ions and L n is the scale length of density profile.For this reason, tokamak edge gyrokinetic codes, e.g.XGC [11,12], are usually limited in the regime where the pedestal width is much greater than the ion gyro-radius.Fully kinetic ion models [13][14][15][16], which retain the ion gyro-motion, avoid these problematic ordering assumptions (although involving more computations in simulations).Low-frequency full kinetics have already been demonstrated in slab geometry [17], successfully benchmarked against gyrokinetics for the slab ITG [18], and extended to toroidal ITG instabilities by Sturdevant et al. [19].In the latter work, a variational integrator was developed to integrate the full orbits of ions in toroidal geometry, which proved to be accurate in capturing both the short-time scale cyclotron motion and long time scale drift motion.However the variational integrator requires solution of a nonlinear equation at each time step, requiring extra computational time.Another limitation of that work is that the equilibrium magnetic field is specified analytically and restricted to circular flux surface configurations.
The present work extends that work in three aspects.First, we implement a simpler full orbit integrator, the Boris integrator [20], in toroidal geometry and demonstrate that the accuracy of this integrator is also sufficient for simulation of ITG instabilities.Second, the equilibrium magnetic configuration is extended to general toroidal configuration specified numerically, enabling simulation of realistic equilibria reconstructed from tokamak experiments.Third, we extend the work in Ref. [19] to the nonlinear regime and investigate the nonlinear saturation of ITG instabilities.To verify the new numerical implementation of the orbit integrator and magnetic configuration, the linear electrostatic ITG frequency and growth rate are compared with those given in Ref. [19] and good agreement is found.The simulation adopts the δ f particle-in-cell (PIC) method.The δ f method is chosen because it is good at resolving perturbations with low amplitude.The PIC method is chosen because it samples the phase-space by the using Monte-Carlo method and the accuracy of the resulting Monte-Carlo integration is superior to the continuum method as dimensionality increases (six dimensional in this work).
The remainder of the paper is organized as follows.Section 2 discusses the fully kinetic ion model for simulating ITG instabilities.The implicit δ f PIC method and the full orbit integrator are briefly discussed in Sec. 3. Section 4 discusses how magnetic configuration and magnetic coordinates are handled in our numerical model.Section 5 gives the linear and nonlinear simulation results of ITG instabilities and the linear benchmark against the results of Ref. [19].A brief summary is given in Sec.
6. Appendix A discusses the electromagnetic model we are working on.

Fully kinetic ion model of ITG instabilities
The fully kinetic ion model for ITG instabilities is described in Ref. [19].The following is a summary of the model.The ion Vlasov equation is written where f i (x, v) is the ion distribution function, x and v are ion position and velocity respectively, E and B are the electric field and magnetic field respectively, q i and m i are the charge and mass of ions respectively.We write f i as an equilibrium part plus a perturbation, i.e., f i = f i0 + δ f i , then Eq. ( 1) is written as where δE and δB are the perturbed part of the electric field and magnetic field, respectively.
To model ITG instabilities, the equilibrium part of ion distribution function f i0 is chosen as [19] where n i0 and T i0 are ion number density and temperature, which depend on a radial variable R r given by where r is the minor radius of magnetic surfaces and b = B 0 /B 0 is the unit vector along the equilibrium magnetic field B 0 .The variable R r is a radial coordinate of the ion guiding-center and thus an approximate constant of motion in the weakly inhomogeneous tokamak magnetic field with equilibrium solution to the kinetic equation (1).Using this form of equilibrium distribution, the kinetic equation ( 2) for the perturbed part of the distribution is written as where κ n i = −n −1 i0 ∂n i0 /∂R r and κ T i = −T −1 i0 ∂T i0 /∂R r are the radial gradients of ion density and temperature, respectively.

Electrostatic limit with adiabatic electrons
In this work, we focus on the electrostatic limit, in which δB is zero and δE = −∇δΦ, where δΦ is the perturbed electric potential.Furthermore, we adopt the simple adiabatic electron model for describing the electron response, in which the perturbed electron density is related to δΦ by where n e0 and T e0 are the equilibrium electron number density and temperature, respectively; e is the elementary charge, . . . is the magnetic surface averaging operator defined by where J is the Jacobian of magnetic coordinates (ψ, θ, φ), ψ is a radial coordinate, θ and φ are poloidal and toroidal angles, respectively.
In the electrostatic limit, Maxwell's equations reduce to Poisson's equation, which further reduces to the quasi-neutrality condition if the space-charge term is neglected.The quasi-neutrality condition is written as where δn i is the perturbed part of the ion number density.Using δn e given by Eq. ( 6) in the above equation, we obtain which serves as our field equation, from which the electric potential δΦ can be solved.We consider modes with n = 0, where n is the toroidal mode number.Then the flux average δΦ is always zero and the field equation ( 9) reduces to an algebraic equation, which can be analytically solved to give We are working on extending the fully kinetic ion model for ITG instabilities to the electromagnetic case with drift-kinetic electrons.This model is discussed in Appendix A.

Implicit δ f particle-in-cell method and full orbit integrator
The ion Vlasov equation ( 5) is solved by using the δ f PIC method [21,22], in which an assembly of markers are loaded in the phase-space according to a distribution function g(x, v).Then the phase space volume occupied by a marker located at (x j , v j ) is given by V ps j = 1/g(x j , v j ).We define the weight of the j th marker by which is the physical particle number carried by δ f i in the phase space volume V ps j .The weight evolution equation is obtained by multiplying both sides of Eq. ( 5) by V ps j and noting that d(V ps j )/dt = 0, yielding where the magnetic perturbation terms have been dropped due to the electrostatic approximation.
An implicit scheme is used to integrate the weight evolution equation.Denoting the right-hand side of the weight evolution equation ( 12) by h(δE, x, v), the implicit scheme we use takes the following form: with We choose the initial guess of δE (n+1) to be equal to δE (n) and then iterate until convergence is achieved.If the iteration is terminated after only two iterations, then this scheme corresponds to the predictor-corrector scheme called Heun's method [23].The field equation ( 9) needs to be solved once in each iteration.
The ion trajectory (x j , v j ) is advanced by a time-centered difference scheme given by for velocity and for position.Here "staggered" time grids are used for v and x: time grids of v are at half-steps while time grids of x are at integer steps.The position at half-steps, x (n+1/2) , which is needed in Eq. ( 13), is approximated by x (n+1/2) = (x (n) + x (n+1) )/2.Further note that the scheme given in Eq (15) is in an implicit form since the unknown v (n+1/2) appears on both sides of the equation.Fortunately, equation ( 15) can be analytically solved in Cartesian basis and its explicit solution is expressed by the Boris algorithm [20].As is discussed in Ref. [24], the above scheme conserves the phase-space volume, which makes it suitable for particle-based methods where phase-space volume conservation is usually implicitly assumed.In this work, cylindrical coordinates are used in integrating the ion orbits.When implementing the Boris scheme in cylindrical coordinates, a local Cartesian coordinate system with basis vectors (e x , e y , e z ) along the local cylindrical basis vectors (e R , e φ , e Z ) at particle location x (n) is set up to perform the velocity integration to obtain v (n+1/2) .Then the particle location is updated in the local Cartesian coordinates by using Eq. ( 16) and then is transformed to the cylindrical coordinates by using the analytic coordinate transformation.After this, the new velocity v (n+1/2) is projected onto the new basis vectors (e R , e φ , e Z ) at particle location x (n+1) .Typical full ion orbits computed by this scheme are compared with the guiding center orbit in Figure 1.This scheme can reproduce correct drift motion even when a large time-step comparable to the gyro-period is used [25].Figure 1 shows examples of orbits computed with large time steps.gyro-period at its initial location (R = 2.1m, Z = 0m, φ = 0).The results show that the full orbits agrees with the guiding-center orbit for the cases with time-step ∆t < T/4.When ∆t is further increased, the computed full orbits deviate from the guiding-center orbit.Further note that the gyro-radius obtained remains nearly the same when the time-step ∆t < T/4.When ∆t is further increased, the gyro-radius becomes larger than the correct value.The magnetic configuration is from EAST tokamak discharge#59954@3.03s.The initial velocity is given by v R = v Z = 1.0 × 10 6 m/s, and v φ = 5 × 10 5 m/s, which corresponds to a kinetic energy of 23 keV.For ∆t = T/16, the orbit is advanced by 23250 time-steps, in which the particle finishes one banana period.
In the linear limit, the fields E and B on the right-hand side of Eq (15) are replaced by the equilibrium fields E 0 and B 0 .In the present work, E 0 = 0 and B 0 is a general toroidal magnetic configuration specified numerically.

Magnetic field specification and field-line-following coordinates
The equilibrium magnetic field is specified numerically by reading and interpolating the output of the equilibrium reconstruction code EFIT [26].This enables us to handle magnetic configurations with arbitrary flux surface shape.
Ion markers are loaded in the field-line-following coordinates (ψ, θ, α).Ion trajectories are integrated in cylindrical coordinates (R, φ, Z) and then linearly interpolated to (ψ, θ, α) coordinates to do the deposition in order to obtain the perturbed ion density at the spatial grids of (ψ, θ, α) coordinates.
At a grid point x k , the perturbed ion density is approximated by where N p is the total number of markers loaded, ∆V s = J (x k )∆ψ∆θ∆α is the volume of the spatial cell, ∆ψ, ∆θ, and ∆α are the grid point spacings in the ψ, θ and α direction respectively, S is the interpolating function defined as with S 1D being the first-order b-spline function given by (then the deposition corresponds to a linear interpolation).In this work, the marker distribution function g is chosen as g = f i0 N p /(V s n i0 ), where V s the spatial volume of the computational box.
After solving the field equation for δΦ, the spatial differential of δΦ is performed in the field-line-following coordinates to determine the perturbed electric field.The deposition and field solving are done in field-line-following coordinates because this coordinate system is efficient for resolving ITG modes, which have k k ⊥ , where k and k ⊥ are the parallel and perpendicular wave-number, respectively.

Simulation results of linear and nonlinear ITG instabilities
In order to benchmark the results against those of Ref. [19], we adopt the DIII-D cyclone base case [29], which is a concentric-circular magnetic configuration.The main parameters used in the benchmarking are summarized in Table .1.
Table 1.DIII-D cyclone base case parameters [29].The safety factor profile is given by q(r) = q 0 + (r − r 0 )q (r 0 ) with q (r 0 ) = ŝq 0 /r 0 , where ŝ is the magnetic shear at r = r 0 (the radial center of the simulation box).In this case R 0 /ρ i = 450.5,where ρ i = v ti /Ω i is the thermal ion gyro-radius at the magnetic axis, v ti = √ T i0 /m i , Ω i = B axis q i /m i is the ion cyclotron angular frequency at the magnetic axis.Deuterium plasma is assumed in our simulation.The radial center of the simulation box is at r = r 0 = 0.24m and the radial width ∆r = 0.11m, which is about 37.5ρ i .The perturbed potential δΦ is set to zero at the radial boundaries.When a marker moves out of the radial boundary, its vertical location is changed from Z to −Z (this is to follow the drift orbit) and its weight is set to zero, where Z is the vertical coordinate of cylindrical coordinates.In the simulation, the perturbed electric potential δΦ is Fourier filtered along the toroidal direction to retain only the n = 29 harmonic, which is further sine transformed along the radial direction and only low radial harmonics are retained.Shown here is the fundamental radial harmonic of δΦ near the low-field-side midplane, which corresponds to k r ρ i = 0.095.δΦ is normalized by T e /e.The frequency and growth rate in this case are ω r /Ω i = 2.388 × 10 −3 and γ/Ω i = 5.8 × 10 −4 , which correspond to the fifth data point in Fig. 4, where the corresponding results from Ref. [19] are ω r /Ω i = 2.423 × 10 −3 and γ/Ω i = 6.0 × 10 −4 .
This is a multi-scale simulation, which includes both the slow-scale ITG instability and the fast-scale wave associated with the ion gyro-motion.In the simulation, we can identify the existence of the ion Bernstein wave (IBW) associated with the ion gyro-motion.The IBW is hidden in the simulation in Fig. 2, the details of which are plotted in Fig. 3 Figure 4 presents the dependence of the linear ITG mode frequency and growth rate on the ion temperature gradient κ T i , which shows that both the frequency and growth rate increase with the temperature gradient drive.Also plotted in Fig. 4 are the fully kinetic results from Ref. [19], which are in good agreement with our new results.Gyrokinetic simulations found that ITG instabilities usually reach peak growth rate near k θ ρ i ≈ 0.3, where k θ is the bi-normal wave-number k θ ≈ nq 0 /r 0 .This trend can also be captured by the fully kinetic ion model.Figure 5 shows the dependence of ITG growth rate and frequency on k θ ρ i given by the fully kinetic model, which shows that the growth rate reaches a peak near k θ ρ i ≈ 0.4.Also plotted in Fig. 5 are the gyrokinetic results from Ref. [30] and [29], which agree with our results within an acceptable discrepancy.The upper horizontal axis shows the corresponding toroidal mode number n.Also plotted are the gyrokinetic results reported in Ref. [30] and [29].Several factors may contribute to the difference between the three results.In our case ρ i = R 0 /450.5, while ρ i = R 0 /500.0 in Ref [30].Safety factor profiles are slightly different: q(r) = q 0 + (r − r 0 )q (r 0 ) in this work, while q(r) = 0.86 − 0.16r/a + 2.52(r/a) 2 in Ref. [30] and q profile variation is neglected in the flux-tube model of Ref. [29].Deuterium plasma is used in our simulation.We furthermore extend the work in Ref [19] to the nonlinear regime and investigate the nonlinear saturation of the ITG instability.Figure 7 plots the nonlinear evolution of the n = 29 ITG instability, which shows that after the linear growth stage, the mode amplitude saturates and then is suppressed.
This suppression is due to the E × B trapping effect [18].We are also studying the nonlinear interaction between ITG modes of multiple toroidal mode numbers and their coupling to the zonal flow, which will be reported in a future publication.

Figure 1 .
Figure 1.Comparison between the full orbits calculated by the Boris scheme with different time step sizes: ∆t = T/16, ∆t = T/8, ∆t = T/4, ∆t = T/2, ∆t = T, and ∆t = 2T, where T is the ion (Deuteron) gyro-period at its initial location (R = 2.1m, Z = 0m, φ = 0).The results show that the full orbits agrees with the guiding-center orbit for the cases with time-step ∆t < T/4.When ∆t is further increased, the computed full orbits deviate from the guiding-center orbit.Further note that the gyro-radius obtained remains nearly the same when the time-step ∆t < T/4.When ∆t is further increased, the gyro-radius becomes larger than the correct value.The magnetic configuration is from EAST tokamak discharge#59954@3.03s.The initial velocity is given by v R = v Z = 1.0 × 10 6 m/s, and v φ = 5 × 10 5 m/s, which corresponds to a kinetic energy of 23 keV.For ∆t = T/16, the orbit is advanced by 23250 time-steps, in which the particle finishes one banana period.

Figure 3 .
Figure 3.Time evolution (a) of the electric potential during tΩ i = [0 : 200] and the corresponding frequency spectrum (b), which shows a clear peak near the ion gyro-frequency.Here f cyc = 2π/Ω i .

Figure 4 .
Figure 4. Dependence of ITG mode frequency and growth rate on the ion temperature gradient κ T i for the DIII-D cyclone base case.The upper horizontal axis shows the ion temperature gradient normalized by the thermal ion gyro-radius.Also plotted are the results from Ref.[19], which are in good agreement with our new results.

Figure 5 .
Figure5.Dependence of ITG mode growth rate and frequency on k θ ρ i for DIII-D cyclone base case.The upper horizontal axis shows the corresponding toroidal mode number n.Also plotted are the gyrokinetic results reported in Ref.[30] and[29].Several factors may contribute to the difference between the three results.In our case ρ i = R 0 /450.5, while ρ i = R 0 /500.0 in Ref[30].Safety factor profiles are slightly different: q(r) = q 0 + (r − r 0 )q (r 0 ) in this work, while q(r) = 0.86 − 0.16r/a + 2.52(r/a) 2 in Ref.[30] and q profile variation is neglected in the flux-tube model of Ref.[29].Deuterium plasma is used in our simulation.

Figure 6 Figure 6 .
Figure6plots the two-dimensional mode structure of the n = 29 ITG instability in the poloidal plane, which shows a clear ballooning structure (i.e., the amplitude on the low-field-side is larger than that on the high-field-side).This is consistent with the physical picture that ITG instabilities are driven by the E × B drift on the tokamak low-field-side while the E × B drift on the high-field side (the "good-curvature" side ) suppresses the instabilities.

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 27 March 2018 doi:10.20944/preprints201803.0221.v1
where L B is the scale length of B 0 .Since arbitrary functions of the constants of motion are solutions to the kinetic equation, the distribution function given by Eq. (3) is approximately an

preprints.org) | NOT PEER-REVIEWED | Posted: 27 March 2018 doi:10.20944/preprints201803.0221.v1 Peer
Although the DIII-D cyclone equilibrium is analytic and circular, the equilibrium is read in as a general equilibrium specified numerically in the G-EQDSK format of EFIT code.No analytic relations particular to this specific configuration is relied on.Preprints (www.-reviewedversion available at plasma 2018, 1, 10; doi:10.3390/plasma1010010