⇦Chinese Physics Letters, 2018, Vol. 35, No. 2, Article code 025201 Particle Trajectory Integrator in Guiding Center Phase Space * Jing-Bo Lin(林竞波)1,2,3, Wen-Lu Zhang(张文禄)2,1,4,3**, Peng-Fei Liu(刘鹏飞)1,2,3, Chao Dong(董超)2,3, Jin-Tao Cao(曹金涛)2,3, Ding Li(李定)2,3,1 Affiliations 1Department of Modern Physics, University of Science and Technology of China, Hefei 230026 2Beijing National Laboratory for Condensed Matter Physics and CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190 3University of Chinese Academy of Sciences, Beijing 100049 4Department of Physics and Astronomy, University of California, Irvine, California 92697, USA Received 30 November 2017 *Supported by the National Natural Science Foundation of China under Grant Nos 11675257, 11675256 and 11705275, the Strategic Priority Research Program of Chinese Academy of Sciences under Grant No XDB16010300, the Key Research Program of Frontier Science of Chinese Academy of Sciences under Grant No QYZDJ-SSW-SYS016, the External Cooperation Program of Chinese Academy of Sciences under Grant No 112111KYSB20160039, and the National Magnetic Confinement Fusion Science Program of China under Grant Nos 2013GB111003 and 2013GB112011.
**Corresponding author. Email: wzhang@iphy.ac.cn Citation Text: Lin J B, Zhang W L, Liu P F, Dong C and Cao J T et al 2018 Chin. Phys. Lett. 35 025201 Abstract A trajectory integrator is developed based on a particle's guiding center Hamiltonian. It is verified by a series of benchmarks, which are in good accordance with theoretical prediction. This integrator can be used as the guiding center trajectory integrator of a particle-in-cell simulation platform, such as the newly developed VirtEx. It can also be used as a stand-alone tool to investigate particle dynamics in a given background field. DOI:10.1088/0256-307X/35/2/025201 PACS:52.65.Cc, 52.55.Fa, 52.65.Tt © 2018 Chinese Physics Society Article Text Calculation of particle motion is the cornerstone of all the particle-in-cell simulation programs. In toroidally confined plasma, the motion of particles is extremely complex due to the toroidal geometry. However, orbital time scale, such as $\tau_{\rm g}$ for gyromotion, $\tau_{\rm d}$ for drift motion and $\tau_{\rm b}$ for bouncing motion of a trapped particle, can be well separated in the presence of strong magnetic field with small inhomogeneity. For phenomena with frequency much lower than gyrofrequency, calculating the particle orbit using conventional fully-kinetic description is not only inefficient but also introducing unnecessary numerical errors. Thus the guiding center description is developed to average out high frequency gyromotion. In the early research by Alfvén,[1] guiding center motion is intuitively separated into the parallel part along the magnetic field line and various perpendicular drift across field lines, $$ \boldsymbol{v}=v_\shortparallel\hat{\boldsymbol{b}}+\frac{c\hat{\boldsymbol{b}}}{eB}\times (mv_\shortparallel^2\hat{\boldsymbol{b}}\cdot\nabla\hat{\boldsymbol{b}}+\mu\nabla B+e\nabla{\it \Phi}).~~ \tag {1} $$ However, this form does not satisfy Liouville's theorem, which is a fatal deficiency for numerical simulation. Boozer et al. introduced a form derived from Liouville's theorem and further developed Hamiltonian theory of guiding center motion[2] using magnetic coordinate. Then White and Chance developed ORBIT/ORBITX[3] code to simulate particle behavior in toroidally confined plasmas. Now efforts have still been made to improving both the theory and code especially for the physics with multiple temporal and spacial scales involved. The relationship between Boozer coordinate, canonical coordinate and flux coordinate is discussed by Meiss and Hazeltine.[4] New coordinates were figured out by White and Zakharov[5] and Wang,[6] which are valid for arbitrary magnetic field while the former Boozer coordinate is valid only in special cases. A series of work on Hamiltonian using canonical coordinate and the concomitant simulations are published by Wang's team.[6-8] A newly developed high-frequency simulation model with fully kinetic ion and gyrokinetic electron provides a clear picture of the relationship between the guiding center's polarization/magnetization and macroscopic flows.[9,10] The concomitant simulation thoroughly discussed the relationship of magnetic moment in configuration space, guiding center space and gyrocenter space.[11,12] Because of the complex toroidal magnetic field structure, it will be convenient to use magnetic coordinates to describe the motion of particles. In the so-called Clebsch representation,[13] the magnetic field lines in 3D space can be written as $$ {\boldsymbol B}=\nabla\alpha\times\nabla\beta, $$ where $\alpha$ and $\beta$ are flux functions with $\nabla\alpha$ and $\nabla\beta$ perpendicular to the magnetic field line. In a realistic toroidal system in fusion plasmas research, magnetic field can be written as $$\begin{align} \boldsymbol B=\,&\nabla(\zeta-q\theta)\times\nabla\psi_{\rm p},~~ \tag {2} \end{align} $$ $$\begin{align} \boldsymbol B=\,&\delta(\psi_{\rm p},\theta)\nabla\psi_{\rm p}+I(\psi_{\rm p})\nabla\theta+g(\psi_{\rm p})\nabla\zeta,~~ \tag {3} \end{align} $$ where the poloidal magnetic flux $\psi_{\rm p}$ is used for labeling magnetic surface, and the exact location on a specific surface is located with poloidal angle $\theta$ and toroidal angle $\zeta$. In this coordinate, the field line is straight with a local helicity $q=d\zeta/d\theta$. Here $g$ is the coil current, $\psi_{\rm p}$ and $I$ are the toroidal flux and plasmas current, respectively. In general, $\delta$ is not zero, but in a device where beta is relatively low and the field has an up-down symmetry, $\delta$ is ignorable.[5] The resulting Jacobian is $$ \mathcal J=[\nabla\psi_{\rm p}\cdot(\nabla\theta \times\nabla\zeta)]^{-1}= \frac{gq+I}{B^2}.~~ \tag {4} $$ To derive the equation of motion for a particle's guiding center, the Hamiltonian could be written as[3] $$ \mathcal{H}=\frac{1}{2}v_\shortparallel^2 + \mu B+{\it \Phi},~~ \tag {5} $$ and the Lagrangian is[4] $$\begin{align} L(\boldsymbol x, \boldsymbol{\dot x},t)=\,&\frac{1}{2}v_\shortparallel(\boldsymbol x, \boldsymbol{\dot x},t)^2+\boldsymbol{\dot x}\cdot \boldsymbol A(\boldsymbol x,t)\\ &-{\it \Phi}(\boldsymbol x,t)-\mu B(\boldsymbol x, t),~~ \tag {6} \end{align} $$ where $\mu=v_\perp^2/(2B)$ is the magnetic moment, and $v_\shortparallel$ is the parallel part of the particle's velocity. For the convenience and conciseness, the time is normalized by the particle's gyrofrequency $\omega_{\rm c}=qB/(mc)$ and the length is normalized by major radius $R_0$ of toroidal geometry hereafter. Then the unit of energy is $m\omega_{\rm c}^2R_0^2$. The magnetic field perturbation in a toroidally confined plasma generally contains both parallel and perpendicular components. The parallel magnetic perturbation, which drives compressional Alfv́en waves and distorts the equilibrium non-resonantly, is not included in our description.[14] The perpendicular magnetic perturbation can thus be written in terms of $\alpha$ as $$ \delta\boldsymbol B=\nabla\times\alpha\boldsymbol B.~~ \tag {7} $$ Thus by choosing $\zeta$ and $\theta$ as the canonical variables, the conjugate canonical momentum reads[4] $$\begin{align} P_\zeta=\,&g(\rho_\shortparallel+\alpha)-\psi_{\rm p},~~ \tag {8} \end{align} $$ $$\begin{align} P_\theta=\,&I(\rho_\shortparallel+\alpha)+\psi,~~ \tag {9} \end{align} $$ where $\rho_\shortparallel=v_\shortparallel /B$ is the normalized parallel velocity. The equation of the particle's motion is $$\begin{align} \dot{P}_\zeta=\,& -\frac{\partial \mathcal{H}}{\partial\zeta}=\rho_\shortparallel B^2\frac{\partial\alpha}{\partial\zeta}-\frac{\partial{\it \Phi}}{\partial\zeta},\\ \dot{P}_\theta =\,& -\frac{\partial \mathcal{H}}{\partial\theta}= -(\mu+\rho_\shortparallel^2 B)\frac{\partial B}{\partial \theta}\\ &+ \rho_\shortparallel B^2\frac{\partial\alpha}{\partial\theta}-\frac{\partial{\it \Phi}}{\partial\theta},\\ \dot{\zeta} =\,& \frac{\partial \mathcal{H}}{\partial P_\zeta}=-(\mu+\rho_\shortparallel^2 B)\frac{\partial B} {\partial\psi_{\rm p}}\frac{I}{D}\\ &+\frac{\rho_\shortparallel^2 B (q+(\rho_\shortparallel+\alpha) I')}{D}\\ &+\rho_\shortparallel^2 B\frac{\partial\alpha}{\partial\psi_{\rm p}}\frac{I}{D}-\frac{I}{D}\frac {\partial{\it \Phi}}{\partial\psi_{\rm p}},\\ \dot{\theta} =\,& \frac{\partial \mathcal{H}}{\partial P_\theta}=(\mu+\rho_\shortparallel^2 B)\frac{\partial B} {\partial\psi_{\rm p}}\frac{g}{D}\\ &+\frac{\rho_\shortparallel^2 B (1-(\rho_\shortparallel+\alpha) g')}{D}\\ &-\rho_\shortparallel^2 B\frac {\partial \alpha}{\partial \psi_{\rm p}}\frac{g}{D}+ \frac{g}{D}\frac{\partial{\it \Phi}}{\partial\psi_{\rm p}}.~~ \tag {10} \end{align} $$ The above equations can be written in a more computation-friendly form by substituting the canonical momentum $P_\theta$ and $P_\zeta$ with Eqs. (8) and (9), $$\begin{alignat}{1} \frac{d\psi_{\rm p}}{dt}=&\frac{ I [ \partial_\zeta{\it \Phi}-\rho_\shortparallel B^2\partial_\zeta\alpha+(\rho^2_\shortparallel B+\mu)\partial_\zeta B]}{D}\\&-\frac{g[ \partial_\theta{\it \Phi}-\rho_\shortparallel B^2 \partial_\theta\alpha+(\rho^2_\shortparallel B+\mu)\partial_\theta B]}{D},~~ \tag {11} \end{alignat} $$ $$\begin{alignat}{1} \frac{d\phi}{dt}=&\frac{\rho_\shortparallel B^2[1-g'(\rho_\shortparallel+\alpha)-g \partial_{\psi_{\rm p}}\alpha]}{D}\\ &-\frac{g[\partial_{\psi_{\rm p}}\alpha{\it \Phi}+(\rho^2_\shortparallel B+\mu)\partial_{\psi_{\rm p}}B] }{D},~~ \tag {12} \end{alignat} $$ $$\begin{alignat}{1} \frac{d\zeta}{dt}=&\frac{1}{D}[-I (\rho^2_\shortparallel B+\mu) \partial_{\psi_{\rm p}}B -I\partial_{\psi_{\rm p}}{\it \Phi}\\ &+\rho_\shortparallel B^2(q+I'(\rho_\shortparallel+\alpha)+I\partial_{\psi_{\rm p}}\alpha)],~~ \tag {13} \end{alignat} $$ $$\begin{alignat}{1} \frac{d\rho_\shortparallel}{dt}=&\frac{1}{D}\{(-1+g'(\rho_\shortparallel+\alpha)+ g\partial_{\psi_{\rm p}}\alpha)[(\rho^2_\shortparallel B\\ &+\mu)\partial_\theta B+\partial_\theta{\it \Phi}]\}\\ &-\frac{1}{D}\{(q+I'(\rho_\shortparallel+\alpha)+I\partial_{\psi_{\rm p}}\alpha)[(\rho^2_\shortparallel B\\ &+\mu)\partial_\zeta B+\partial_\zeta {\it \Phi}]\}\\ &-\frac{1}{D}\{(g\partial_\theta\alpha-I\partial_\zeta\alpha)[\partial_{\psi_{\rm p}}B+\partial_{\psi_{\rm p}}{\it \Phi}]\}\\ &-\partial_t\alpha.~~ \tag {14} \end{alignat} $$ These are the equations of motion evaluated in the code. In a toroidal geometry, particles may produce two kinds of trajectories by the shape of trajectories projected on the cross section. A passing particle's trajectory projection to the poloidal plane is almost a circle while a trapped particle's trajectory projection is like a banana. For a passing particle, the motion of the guiding center can be separated into two parts, parallel motion $v_\shortparallel$ along the magnetic field line and drift motion,[15] $$ \boldsymbol v_{\rm d}=\frac{W+W_\shortparallel}{B_\zeta R}\hat{\boldsymbol{n}}\times\hat{\boldsymbol{b}},~~ \tag {15} $$ where $B_\zeta=\boldsymbol{B} \cdot \nabla\zeta / |\nabla\zeta|$ is the toroidal magnetic field, $R$ is the major radius, $W_\shortparallel=v_\shortparallel^2/2$ and $W_\perp=\mu B$ are the particle's kinetic energy parallel and perpendicular to magnetic field, respectively, $W=W_\shortparallel+W_\perp$ is the total kinetic energy, $\hat{\boldsymbol{n}}$ is the norm vector of the field line curvature, and $\hat{\boldsymbol{b}}=\boldsymbol B/B$ is the unit vector in magnetic field direction. Thus the guiding center can be described by $$\begin{align} \frac{dR}{dt}\approx\,&\frac{v_\shortparallel}{qR_0}Z,~~ \tag {16} \end{align} $$ $$\begin{align} \frac{dZ}{dt}\approx \,&- \frac{v_\shortparallel}{qR_0}(R-R_0) + v_{\rm d},~~ \tag {17} \end{align} $$ where $Z$ is one of the cylindrical axes perpendicular to the toroidal plane, and $R_0$ is the major radius of magnetic axis. The particle's guiding center trajectory $(R-R_0-v_{\rm d}qR_0/v_\shortparallel)^2+Z^2\simeq {\rm const}$ is almost a circle with center shifted from the original magnetic axes. The shift distance is $$ d=v_{\rm d}qR_0/v_\shortparallel.~~ \tag {18} $$ For a trapped particle, since the inside magnetic field strength is stronger than the outside, the particle is trapped in a distorted magnetic mirror. With the mirror force expressed as $-\mu({dB}/{ds})$, the equation of motion is $$ \frac{d^2s}{dt^2}=-\mu\frac{dB}{ds},~~ \tag {19} $$ where $s$ is the distance along a magnetic field line, and $dB/ds$ can be evaluated using the geometric relation $B_\phi ds=Brd\theta$. The solution reads $$ \theta=\Big(\frac{2W_\shortparallel R_0}{W_\perp r}\Big)^{1/2} \sin\Big[\Big(\frac{rW_\perp}{R_0^2q^2}\Big)^{1/2} t\Big].~~ \tag {20} $$ Thus the trapped particle is bouncing in the $\theta$ direction with amplitude $\theta_t=[{2W_\shortparallel R_0}/({W_\perp r})]^{1/2}$ at frequency $\omega_t=[{rW_\perp}/({R_0^2q^2})]^{1/2}$. The derivative of $\theta$ with respect to time is $$ \frac{d\theta}{dt}=\omega_t\cos(\omega_t t).~~ \tag {21} $$ The motion on $r$ direction is dominated by the curvature drift $$ \frac{dr}{dt}=v_{\rm d}\sin\theta\approx\frac{W_\perp}{R_0B_0}\theta.~~ \tag {22} $$ The parameter equations of the particle's trajectory can be obtained using Eqs. (21) and (22), $$\begin{align} &\frac{d\theta}{dr}=\frac{{R_0B_0}\theta_t\omega_t(1-\theta^2/\theta_t^2)^{1/2}} {{W_\perp}\theta},~~ \tag {23} \end{align} $$ $$\begin{align} &\Big(\frac{\theta_tW_\perp}{R_0B_0}\Big)^2(1-\theta^2/\theta_t^2) =\omega_t^2(r-r_0)^2,~~ \tag {24} \end{align} $$ where $r_0$ is a constant as the result of the integration, which represents the center of the 'banana'. The width of the banana is $$ \Delta r=2{\theta_tW_\perp}/({\omega_t R_0B_0}).~~ \tag {25} $$ Based on Eqs. (11)-(14), a numerical code was developed to calculate the particle's guiding center trajectory in toroidal geometry. In the code, background field properties are interpolated from grid point using second order interpolation. The field information of a certain equilibrium can be read from an external file. Moreover, the program itself can also generate a simple equilibrium for a high aspect ratio tokamak scenario. Diagnostic routines are also added for tracing and checking the energy conservation. To ensure the fidelity of this pusher code, we set some simple scenarios to benchmark this program. We first benchmarked individual drift motions. The theoretical drift speeds of the particle's guiding center are[16] $v_{R}=2W_\shortparallel(\boldsymbol{R}_{\rm c}\times\boldsymbol{B})/(q R_{\rm c}^2 B^2)$ for curvature drift, where $\boldsymbol{R}_{\rm c}$ is the curvature radius vector, $v_{\nabla B}=\pm \mu (\boldsymbol{B}\times\nabla B)/(qB^2)$ for $\nabla B$ drift, and $v_{E}=\boldsymbol{E}\times \boldsymbol{B}/B^2$ for $\boldsymbol{E}\times\boldsymbol{B}$ drift. Figure 1 shows the example $\boldsymbol{E}\times\boldsymbol{B}$ drift trajectory of the guiding center. Figures 2–4 show the calculation results obtained in several benchmarks including $\boldsymbol{E}\times \boldsymbol{B}$ drift, $\nabla B$ drift and curvature drift, respectively. In $\boldsymbol{E}\times\boldsymbol{B}$ drift test, both magnetic field $\boldsymbol{B}$ and electric field $\boldsymbol{E}$ were constant and homogeneous. In $\nabla B$ drift test, the background magnetic field was set to $\boldsymbol{B}=(1/R) \hat{\boldsymbol\zeta}$, where $\hat{\boldsymbol\zeta}=\nabla\zeta/|\nabla\zeta|$ is the unit vector on the $\zeta$ direction. The particle's parallel velocity $v_\shortparallel$ in $\nabla B$ drift test was set to 0 to elude the interference of curvature drift. In curvature drift, the background magnetic field was also set to $\boldsymbol{B}=1/R \hat{\boldsymbol{\zeta}}$, while the perpendicular kinetic energy was set to $W_\perp=0.05W$ to minimize the effect of $\nabla B$ drift. In these three cases, the code's results match the theoretical curves well. It is notable that in the curvature drift case, the $\nabla B$ drift effect was also included in the theoretical curve.
cpl-35-2-025201-fig1.png
Fig. 1. Guiding center trajectory of $\boldsymbol{E}\times \boldsymbol{B}$ drift. The background electric field is set to $1\times 10^{-6}$ in normalized unit, and the simulation time duration is 50 ion gyroperiod.
cpl-35-2-025201-fig2.png
Fig. 2. Guiding center drift distance benchmark for $\boldsymbol{E}\times\boldsymbol{B}$ drift.
As an integrated test, particle motion was benchmarked in the tokamak geometry. Here for simplicity, a high aspect ratio tokamak with low $\beta$ limit was chosen. In this scenario, we may choose the reverse aspect ratio $a/R$ as a smallness factor, then the magnetic field can be simplified to[3] $$ \boldsymbol B=[r/q(r)]\nabla\phi\times\nabla r+R_0\nabla\phi,~~ \tag {26} $$ where $\phi$ is a geometrical toroidal angle, $\theta$ is a poloidal angle, and $r$ is the minor radius. The Shafranov shift can be ignored for a concentric cross section. In the following benchmark, $R_0=1$, $B(R_0)=1$ and $\psi_{\rm pw}=0.02895$, where $\psi_{\rm pw}$ is the poloidal magnetic flux measured on the conductive metal wall. The safety factor is shown in Fig. 5.
cpl-35-2-025201-fig3.png
Fig. 3. Guiding center drift distance benchmark for $\nabla{B}$ drift.
cpl-35-2-025201-fig4.png
Fig. 4. Guiding center drift distance benchmark for curvature drift.
cpl-35-2-025201-fig5.png
Fig. 5. Safety factor used in benchmark, and $q=1.28+0.80\psi_{\rm p}+1.00\psi_{\rm p}^2$.
Two kinds of particle were benchmarked: trapped and passing particles, respectively. In both the cases, the parameters of particle are $m=1$ and $q=1$ in the normalized unit. The initial position of particle was $\psi=0.30\psi_w$, $\theta=0$ and $\zeta=0$. For the trapped particle, the ratio between perpendicular kinetic energy and total energy was 0.98. The banana width of the trapped particle was measured in various cases, in which the particle's normalized energy was set from $1\times 10^{-6}$ to $4 \times 10^{-6}$. The trapped particle's trajectory projection on a poloidal cross section, particle's trajectory on $\zeta$ direction, relative energy error result and comparison with theoretical estimation Eq. (25) are shown in Fig. 6. The relative energy error plot shows that, for several bounce periods, the relative error was of $10^{-6}$ order, which is acceptable for most applications.
cpl-35-2-025201-fig6.png
Fig. 6. (a) The trajectory projected on cross section. The cross symbol + indicates the center of the banana orbit trajectory. (b) Relative energy error, where $W_0$ is particle's normalized initial energy, and $W$ is the normalized particle energy. (c) Toroidal motion of a trapping particle. (d) Benchmarks with various particle energies.
cpl-35-2-025201-fig7.png
Fig. 7. (a) The trajectory projected on cross section. The cross symbol + indicates the center of the guiding center drift trajectory. (b) Relative energy error. (c) Toroidal motion of a trapping particle. (d) Benchmarks with various particle energies.
For the passing particle, the ratio between perpendicular kinetic energy and total energy was 0.05. The distance between the drift orbit center and the magnetic axis was compared with theoretical expression (18). The simulation result and comparison with theoretical estimation are shown in Fig. 7. The relative error of energy showed the error level was also at $10^{-6}$ level. In both the trapped and passing particles, the subtle discrepancy between theoretical calculation and the simulation result was due to the approximation made in theoretical derivation. In summary, a particle trajectory integrator calculating the motion of the particle's guiding center in toroidal geometry has been reported. We first introduce the magnetic field coordinate, which is widely used in fusion plasmas. By using magnetic coordinate, the guiding center Hamiltonian, canonical variables and momentum are derived. Based on the guiding center motion equation, we develop a particle trajectory integrator to calculate the particle's behavior inside toroidally confined plasmas. The benchmark of this code shows a good agreement in high-aspect ratio magnetic field with two typical types of particles, whose behavior is already well known. This code can be used not only in self-consistent simulations in plasma simulation platform, such as VirtEx,[17] but also as a stand-alone tool for investigating various particle transport and acceleration problems in fusion and space plasmas.[18,19] This research used resources of the National Supercomputer Center in Tianjin, the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory and the National Energy Research Scientific Computing Center. References Guiding center drift equationsHamiltonian guiding center drift orbit calculation for plasmas of arbitrary cross sectionCanonical coordinates for guiding center particlesHamiltonian guiding center equations in toroidal magnetic configurationsCanonical Hamiltonian theory of the guiding-center motion in an axisymmetric torus, with the different time scales well separatedExplicit Runge–Kutta integrator with Hamiltonian correction for long-time simulations of guiding-center orbit in tokamak configurationsLinear gyrokinetic theory and computation of the gyrocenter motion based on the exact canonical variables for axisymmetric tokamaksA closed high-frequency Vlasov–Maxwell simulation model in toroidal geometryPushforward transformation of gyrokinetic moments under electromagnetic fluctuationsParticle simulation of radio frequency waves with fully-kinetic ions and gyrokinetic electronsHamiltonian theory of guiding-center motionDestruction of magnetic surfaces by magnetic field irregularitiesVerification of linear resistive tearing instability with gyrokinetic particle code VirtExElectron Acceleration by Small-Amplitude Solitary Kinetic Alfvén Wave in a Low-Beta PlasmaElectron Acceleration by a Finite-Amplitude Solitary Kinetic Alfvén Wave
[1]Alfvén H 1940 Ark. För Matematik Astronomi Och Fys. 27A 1
[2] Boozer A H 1980 Phys. Fluids 23 904
[3] White R B and Chance M S 1984 Phys. Fluids 27 2455
[4] Meiss J D and Hazeltine R D 1990 Phys. Fluids B: Plasma Phys. 2 2563
[5] White R B and Zakharov L E 2003 Phys. Plasmas 10 573
[6] Wang S 2006 Phys. Plasmas 13 052506
[7] Xiao X T and Wang S J 2008 Phys. Plasmas 15 122511
[8] Xu Y F, Xiao X T and Wang S J 2011 Phys. Plasmas 18 042505
[9] Liu P F, Zhang W L, Dong C, Lin J B, Lin Z H and Cao J T 2017 Nucl. Fusion 57 126011
[10] Liu P F, Zhang W L, Dong C, Lin J B, Lin Z H, Cao J T and Li D 2017 Phys. Plasmas 24 112114
[11] Lin J B, Zhang W L, Liu P F, Lin Z H, Dong C, Cao J T and Li D 2018 Nucl. Fusion 58 016024
[12]Lin J B, Zhang W L, Liu P F, Cao J T and Li D 2018 Commun. Comput. Phys. (submitted)
[13] Cary J and Brizard A 2009 Rev. Mod. Phys. 81 693
[14] Rosenbluth M N, Sagdeev R Z, Taylor J B and Zaslavski G M 1966 Nucl. Fusion 6 297
[15]Wesson J 2011 Tokamaks (Oxford: Oxford University Press)
[16]Chen F F 1984 Introduction to Plasma Physics and Controlled Fusion (London: Springer)
[17] Feng H Y, Zhang W L, Dong C, Cao J T and Li D 2017 Phys. Plasmas 24 102125
[18] Wang D Y, Song Q W and Yang L 2008 Chin. Phys. Lett. 25 794
[19] Wang D Y 2009 Chin. Phys. Lett. 26 019601