Global gyrokinetic particle simulation of resistive tearing modes has been developed and verified in the gyrokinetic toroidal code (GTC). GTC linear simulations in the fluid limit of the kink-tearing and resistive tearing modes in the cylindrical geometry agree well with the resistive magnetohydrodynamic eigenvalue and initial value codes. Ion kinetic effects are found to reduce the radial width of the tearing modes. GTC simulations of the resistive tearing modes in the toroidal geometry find that the toroidicity reduces the growth rates.

Since the early 1960s, macroscopic tearing modes have drawn widespread attention from plasma physicists in both space and laboratory plasmas. In space plasmas, tearing mode is considered to be the underlying mechanism for the solar flare in the sun1 and the substorm in the earth magnetosphere.2,3 In laboratory plasmas, especially for magnetic fusion plasmas, macroscopic tearing mode may cause the sawtooth oscillations and induce disruptions in tokamak discharge.4 Tearing mode is driven by plasma current and will release the free energy of non-uniform magnetic field. In fusion experiments, one can usually control the equilibrium current profile to avoid the q = 1 rational surface, where the parallel wave number k = 0 for (m = 1, n = 1) mode; here, q is the safety factor, m is the poloidal mode number, and n is the toroidal mode number. In this way, the most dangerous (1, 1) tearing mode can be avoided. However, other macroscopic tearing modes, for example, (2, 1) and (3, 2) modes, could remain unstable. Even if the classical tearing modes are stable, their neoclassical version, neoclassical tearing modes, could still be unstable and pose a severe threat to the plasma confinement. So the control of the tearing modes is crucial to the future fusion experiment in ITER5 and need to be carefully investigated.

The tearing mode has been studied in fusion plasmas for several decades.6,7 However, important physics such as kinetic effects on the tearing mode remains as an unsolved problem. The reason is that tearing mode is characterized with both multi-spatial scale, which spans from electron skin depth De to tokamak size, and multi-time scale, which varies from plasma Alfvén time τA to collisional dissipation time 1/νei. Because of the multiscale nature, it is difficult to make accurate theoretical analysis of the tearing mode. Theorists typically separate the problem into two regions, an ideal magnetohydrodynamic (MHD) region and an inner region, where the non-ideal effects, such as plasma resistivity and electron inertial,8,9 play the role as the releasing channel of free magnetic energy. By matching the inner and outer solutions, the method can provide an estimate of the dispersion relation of the tearing mode.6 Some of the simulation work of the tearing mode have been performed on large parallel simulation codes, such as M3D10 and NIMROD.11 In these codes, the background plasma is modeled as a resistive MHD, and the energetic particles are simulated as guiding centers. Indeed, using these codes, one can investigate the fluid behavior of the tearing mode and study the kinetic effects of energetic particles on the tearing mode.

In the high temperature tokamak discharge, kinetic effects, such as orbital effects of both thermal and energetic particles, will become important and affect the behavior of the tearing mode. The time and spatial scales of particle characteristic motion will give rise to more complexity in the tearing mode study.12,13 To accurately predict the behavior of the tearing mode in fusion plasmas, one needs to start from the more realistic first-principles physics model in the tokamak geometry, which is the motivation of this work. We utilize the Gyrokinetic Toroidal Code (GTC),14,15 which has been extensively applied to study neoclassical and turbulent transport,16,17 energetic particles,18 and Alfvén eigenmodes.19,20 The implementation of the equilibrium current21 enables simulation of current-driven instabilities in the toroidal geometry with kinetic effects, such as the kink mode22 and the tearing mode in this paper.

To efficiently treat the electron dynamics in the simulation of the ion scale turbulence, a fluid-kinetic hybrid electron model has been implemented in GTC,23 by assuming ω/kvte 1; here, ω is the mode frequency of interest, k is the parallel wave number, and vte is the thermal speed of electrons. It can efficiently simulate the low frequency MHD modes in the lowest order massless fluid model, and calculate the wave-electron interaction in the higher orders by solving the drift kinetic equation. In order to recover the tearing mode, we need to resolve the singularity near the mode rational surface via some non-ideal MHD effects, such as the resistivity or finite electron inertia. In this paper, we will focus on the classical tearing mode excited by the resistivity. We first introduce a resistive electron fluid model that can recover the classical tearing mode in Sec. II. We then verify GTC simulation results of the classical tearing mode using an eigenvalue code in cylinder in Sec. III. Finally, we incorporate ion kinetic effects on the tearing mode and extend the simulation of the tearing mode to the realistic tokamak geometry in Sec. IV. GTC simulation of collisionless tearing mode driven by the electron inertia has also been verified and will be reported in a follow-up paper.

In the resistive tearing mode, the width of the perturbed current layer increases with larger resistivity. In this work, we assume that the resistivity is sufficiently large such that the width of the current layer is much bigger than the electron gyroradius. The electron dynamics can, thus, be described by the drift kinetic equation. In guiding center coordinates (R, μ, v), time evolution of electron distribution function fe follows:

ddtfe(R,μ,v||,t)=(t+Ṙ+v̇||v||)fe=(tfe)collision,
(1)

where R is the position of electron guiding center, μ and v are the magnetic moment and velocity along the magnetic field of electron. Here, we will use the Krook collisional operator (tfe)collision=νei(fefe0) in the following derivation, where fe0 is the equilibrium distribution function of electron. The drift velocity of electron guiding center, which includes the parallel, electric, magnetic gradient, and curvature drifts, reads

dRdt=v||BB0+cb0×ϕB0+v||2Ωe×b0+μmeΩeb0×B0.
(2)

The parallel acceleration due to the mirror force and the electric field is

dv||dt=1meB*B0(μB0+qeϕ)qemecA||t,
(3)

where

B*=B0+B0Ωe×b0+δB.
(4)

Here, R, v, μ, me, qe, and Ωe denote guiding center position, the parallel velocity, the magnetic moment, the mass, the electric charge, and cyclotron frequency of the electron, and B0, δB, ϕ, and A denote the equilibrium magnetic field, the perturbed magnetic field, the electrostatic potential, and parallel vector potential, respectively. b0=B0B0 is the unit vector of equilibrium magnetic field, c and t denote the light speed and the time. Assuming there is no equilibrium electric field, ϕ and A can be replaced by their perturbed part δϕ and δA, respectively, then δB can be written as δB=×δAb0 approximately. Assuming a shifted Maxwellian for the equilibrium electron distribution function that satisfies the 0th order electron drift kinetic equation, we can calculate the moments of the drift kinetic equation order by order, and derive the perturbed fluid continuity equation of electron21 

tδne+B0(ne0δu||eB0)+B0δvE(ne0B0)ne0(δv*e+δvE)B0B0+δB(ne0u||e0B0)+c×B0B02(δpee+ne0δϕ)+{δB(ne0δu||eB0)+B0δvE(δneB0)+c×B0B02δneδϕ}NL=0
(5)

and the parallel momentum equation

ne0metδu||e+ne0meu||e0δu||eb0=ne0e(||δϕ1cδA||t)δBB0pe0||δpene0meνeiδu||e.
(6)

Here, the total electron density ne, parallel flow ue, and pressure pe are the sum of their equilibrium and perturbed parts, i.e., ne=ne0+δne, ue=ue0+δue, and pe=pe0+δpe. In the ion frame, the current is carried by parallel electron flow. We ignore the ion flow in this work, i.e., the simulation is performed in the plasma frame ignoring the flow shear effects. Then, the equilibrium Ampère law can be written approximately as

ne0qeu||e0J||0=4πc(×B0)||.
(7)

To close the system of the fluid equation, we use the isothermal electron model, i.e., Te = constant, pe=neTe, then

δpe=δneTe+δr(ne0Te),
(8)

where δr is the displacement of the fluid elements. One can also apply a more rigorous energy equation for electron, but needs to calculate the higher order moment from the electron drift kinetic equation.15 Since we are focusing on resistive tearing mode in this work, we will neglect the electron inertial term, which will be discussed in the follow-up paper on collisionless tearing mode. Then, we can use the massless electron momentum equation to evolve the parallel vector potential as

tδA||=cb0δϕ+cne0eb0δpe+cne0eδBB0pe0+meceνeiδu||e.
(9)

The electrostatic potential is calculated from the gyrokinetic Poisson's equation24 

4πZi2Ti(δϕδϕ̃)=4π(Ziδnieδne).
(10)

Here, the left-hand-size of Eq. (10) represents the ion polarization density and δϕ̃ is the second gyrophase-averaged potential as defined in Ref. 24.

The electron flow is calculated from the parallel Ampère's Law

δu||e=c4πene02δA||+Zini0ene0δu||i.
(11)

In these equations, the ion guiding center density ni and current ui can be calculated from the standard gyrokinetic model.25 

The fluid electrons (5), (8), (9) and the gyrokinetic ions are coupled through the gyrokinetic Poisson's equation (10) and Ampère's law (11). These equations form a closed system, which can simulate the low frequency electromagnetic MHD waves and the interactions between these modes and particles.

As the first step to investigate the tearing mode, we need to recover the MHD behavior of the tearing mode via the resistive fluid part of this model and verify the new capability of GTC. Then, we will study kinetic ion effects on tearing mode with the gyrokinetic ions. We will use magnetic flux coordinates (Ψ, θ, ζ) convenient for tokamak, where Ψ is the poloidal magnetic flux, θ and ζ are, respectively, poloidal and toroidal angles. Using this coordinate system, the equilibrium magnetic field26 can be represented as either covariant form

B0=δΨ+Iθ+gζ
(12)

or contravariant form

B0=qΨ×θΨ×ζ,
(13)

and the Jacobian is

J1=Ψ(θ×ζ)=B02gq+I.
(14)

The low-mode-number tearing modes, such as, (1, 1) and (2, 1) tearing modes are very dangerous macroscopic MHD instabilities in tokamak discharge. We will focus on these two modes in the following discussions. To verify the capability of GTC for tearing mode simulation, we have developed a simple 1D (in minor radius direction) initial value code and an eigenvalue code using resistive MHD in the cylindrical geometry.

Applying an initial perturbation: δA||=δÂ||(Ψ)cos(mθnζ), and calculating the m and n harmonics of perturbed quantities, we can get a simplified one-dimensional model.21 Expending the ϕ̃ with the modified Bessel function and neglecting the gyrokinetic ion contribution, the gyrokinetic Poisson's equation (10) in the long wave length limit approximation ρi221 becomes

ρs2λDe22δϕ=4πδneqe.
(15)

Here, ρs is the ion gyroradius defined using ion sound wave speed and λDe is the electron Debye length.

For simplicity, we will consider the uniform equilibrium pressure and keep the equilibrium current driven term for tearing mode. As electrons E×B drift cancels with that of ions, the electron continuity equation can be simplified to

δnet+B0(ne0δu||eB0)+δB(ne0u||e0B0)=0.
(16)

Together with the massless electron parallel momentum equation

δA||t=cb0δϕ+cne0eb0δpe+meceνeiδu||e
(17)

and parallel Ampère's law

δu||e=c4πene02δA||,
(18)

where 2=1/r((r/r)m2/r2), ||=(nm/q)BT/RB0, and δB=×δAb0. Using this 1D model, an initial value code and an eigenvalue code have been developed. Before the simulation of tearing mode, we reproduced (1, 1) kink mode in the cylindrical geometry, which has already been studied via gyrokinetic simulation.22,27,28

The (1, 1) tearing mode is also called the kink-tearing mode because of the strong coupling between the kink mode and the tearing mode. The resistivity can increases the growth rate of the (1, 1) ideal kink mode. We use a q profile, q=0.8+3.2(r/R0)2, where r is the minor radius. The other simulation parameters are: inverse aspect ratio ε = a/R0 = 0.5, major radius R0 = 100 cm, magnetic field B0 = 10 000G, equilibrium electron density on magnetic axis ne0 = 1014/cm3, and equilibrium electron beta βe = pe0/(8πB02) = 0.004. Note that the position of q = 1 mode rational surface is r = 0.25R0. For ideal MHD, in which δA||/t=cb0δϕ, we get the ideal kink mode structures in Fig. 1.

FIG. 1.

Radial mode structures of (1, 1) kink mode from both GTC simulation and eigenvalue calculation in the cylindrical geometry. The position of the mode rational surface q = 1 is at r/R0 = 0.25.

FIG. 1.

Radial mode structures of (1, 1) kink mode from both GTC simulation and eigenvalue calculation in the cylindrical geometry. The position of the mode rational surface q = 1 is at r/R0 = 0.25.

Close modal

In Fig. 1, the red and blue curves are the radial mode structures of δA and δϕ, respectively. The radius mode structures are measured at θ = 0 and each of the curves is normalized to their maximum values, since the absolute values of the mode amplitudes are not important in the linear simulations. The growth rate of eigenvalue method γeig and growth rate of GTC γgtc equal to 0.043ωA and 0.038ωA, respectively. Because of the ideal MHD approximation, the mode structure of δA is always 0 and δϕ is steep on the mode rational surface with q = 1 (i.e., r/R0 = 0.25), which are consistent with the typical mode structures of (1, 1) kink mode.22,27

When the equilibrium current profile satisfies the criteria Δ > 0, the classical tearing mode is unstable. Here, Δ is defined as the logarithmic jump of the perturbed magnetic flux in the outer ideal MHD region.6,11Δ=[δψ(rs+δr)/rδψ(rs+δr)/r]/δψ(rs); here, rs is the radius position of rational surface, and δr is the width of the inner resistive region. The (1, 1) ideal kink mode is a typical ideal MHD mode driven by equilibrium current. When the resistivity is considered in the electron parallel momentum equation, Eq. (9), (1, 1) kink mode will appear as the kink-tearing mode.

Due to the finite resistivity, there will be finite parallel electric field on the mode rational surface where k(r)=0. In contrast to the ideal MHD kink mode, the parallel inductive vector potential δA is not zero on the mode rational surface, which can be seen in Fig. 2 for the νei = 0.02Ωp kink-tearing case. The finite δA induces the finite δBr on the rational surface, which reconnects the field line and tears the magnetic surface. Because the finite resistivity will resolve the singularity of ideal MHD, the mode structure of the tearing mode will be broader near the mode rational surface, which can also be seen in Fig. 2, compared with the mode structure of the ideal kink mode. The half width of the perturbed current for the tearing mode is about 0.1R0, which is much larger than the thermal ion gyroradius 0.0022R0 for the above parameters. In Fig. 3, the left and right panels are the (1, 1) kink-tearing mode structures of δA|| and δϕ. They are typical (1, 1) kink-tearing mode structure and show poloidal symmetry in cylinder simulation. For the specific νei = 0.02Ωp case, the growth rate of tearing mode from eigen value, GTC and 1D initial value are γeig = 0.074 ωA, γgtc = 0.066ωA, and γinit = 0.074ωA. Although the growth rate of (1, 1) kink-tearing mode is larger than the pure kink mode, the growth rate of the (1, 1) tearing mode is still dominated by the kink mode.

FIG. 2.

Comparison of the radial mode structures between (1, 1) ideal kink and kink-tearing mode from GTC simulation (-gtc represents GTC results), 1D eigenvalue (-eigen represents eigenvalue results), and initial value (-1D represents initial value results) calculations in the cylindrical geometry with the mode rational surface of q = 1 at r = 0.25R0. Panel (a) is for parallel vector potential and panel (b) is for electron current perturbation.

FIG. 2.

Comparison of the radial mode structures between (1, 1) ideal kink and kink-tearing mode from GTC simulation (-gtc represents GTC results), 1D eigenvalue (-eigen represents eigenvalue results), and initial value (-1D represents initial value results) calculations in the cylindrical geometry with the mode rational surface of q = 1 at r = 0.25R0. Panel (a) is for parallel vector potential and panel (b) is for electron current perturbation.

Close modal
FIG. 3.

Poloidal contour plots of δA|| (panel (a)) and δϕ (panel (b)) mode structure of (1, 1) kink-tearing mode from GTC simulation in the cylindrical geometry.

FIG. 3.

Poloidal contour plots of δA|| (panel (a)) and δϕ (panel (b)) mode structure of (1, 1) kink-tearing mode from GTC simulation in the cylindrical geometry.

Close modal

Furth et al.6 have predicted theoretically that the growth rate γ depends on η3/5 for a pure tearing mode. To verify the GTC capability of the tearing mode simulation, one needs to recover this prediction for the pure tearing mode. Since the (2, 1) kink mode is always stable when η = 0, we investigate the (2, 1) tearing mode behavior in the cylinder with q=1.61550.327r/R0+4.232(r/R0)2, ε = 0.58, and all other simulation parameters same as those of the above (1, 1) mode simulation. For a specific νei = 0.001Ωp, the (2, 1) tearing mode structures are shown in Fig. 4.

FIG. 4.

Comparison of the mode structures of (2, 1) tearing mode from GTC simulation and 1D eigenvalue calculation in the cylindrical geometry. The position of mode rational surface q = 2 is at r/R0 = 0.343.

FIG. 4.

Comparison of the mode structures of (2, 1) tearing mode from GTC simulation and 1D eigenvalue calculation in the cylindrical geometry. The position of mode rational surface q = 2 is at r/R0 = 0.343.

Close modal

In Fig. 4, the solid lines and dotted lines are the radial mode structures from eigenvalue calculation and GTC simulation, respectively. All these mode structures are typical (2, 1) tearing mode structure7,29,30 and the growth rate from GTC and eigenvalue method are γeig = 0.0050ωA and γgtc = 0.0043ωA, respectively. In GTC simulation, the poloidal mode structures of both δA|| and δϕ are typical (2, 1) tearing mode structure as shown in Fig. 5.

FIG. 5.

Poloidal contour plots of δA|| (panel (a)) and δϕ (panel (b)) mode structure of (2, 1) tearing mode from GTC simulation in the cylindrical geometry.

FIG. 5.

Poloidal contour plots of δA|| (panel (a)) and δϕ (panel (b)) mode structure of (2, 1) tearing mode from GTC simulation in the cylindrical geometry.

Close modal

In Figs. 4 and 5, the GTC simulations have accurately recovered the (2, 1) tearing mode structures. The growth rate dependence of (2, 1) tearing mode on resistivity agrees well with the theoretical prediction as in Fig. 6.

FIG. 6.

The dependence of the (2, 1) tearing mode growth rate γ on the resistivity η.

FIG. 6.

The dependence of the (2, 1) tearing mode growth rate γ on the resistivity η.

Close modal

It can be seen that the simulation result is consistent with the theoretical prediction when the resistivity is small. Small resistivity indicates large time scale separation between the tearing mode growth time 1/γ and the dissipation time 1/νei. With this time scale separation, theoretical prediction works properly. However, for the excessively lager resistivity, the modes will also gradually be damped by the collisional dissipation and the growth rate scaling deviates from the η3/5 dependence.

One of the important issues for tearing modes, which has not been adequately addressed in fusion plasmas, is the thermal ions kinetic effects. Calculating the integral of perturbed distribution function of gyrokinetic ions, we can get the charge and current density of thermal ions. Then, substitute them into the Poisson's equation and parallel Ampère's Law, we can study the thermal ion kinetic effects on the tearing mode. Using the same parameters as those in the fluid simulations with Ti = Te = 500 ev and βe = 0.004, but treating thermal ions as gyrokinetic or drift kinetic, we simulate the (1, 1) kink-tearing mode using νei = 0.02Ωp and (2, 1) tearing mode using νei = 0.001Ωp. The mode structures and growth rates are shown in Fig. 7 and Table I for the (1, 1) kink-tearing mode and in Fig. 8 and Table II for the (2, 1) tearing mode.

FIG. 7.

The mode structures of (1, 1) kink-tearing mode in the cylindrical geometry from GTC simulations with different thermal ions: fluid ion (iMHD), drift kinetic ions (DKI), and gyrokinetic ions (GKI).

FIG. 7.

The mode structures of (1, 1) kink-tearing mode in the cylindrical geometry from GTC simulations with different thermal ions: fluid ion (iMHD), drift kinetic ions (DKI), and gyrokinetic ions (GKI).

Close modal
TABLE I.

Growth rate of (1, 1) kink-tearing mode with different thermal ions from GTC simulations in the cylindrical geometry.

Resistive MHDFluid ionDKIGKI
0.060ωA 0.061ωA 0.0467ωA 0.0467ωA 
Resistive MHDFluid ionDKIGKI
0.060ωA 0.061ωA 0.0467ωA 0.0467ωA 
FIG. 8.

The mode structures of (2, 1) tearing mode in the cylindrical geometry from GTC simulation with fluid, drift kinetic ions, and gyrokinetic ions.

FIG. 8.

The mode structures of (2, 1) tearing mode in the cylindrical geometry from GTC simulation with fluid, drift kinetic ions, and gyrokinetic ions.

Close modal
TABLE II.

Growth rate of (2, 1) tearing mode with different thermal ions from GTC simulations in the cylindrical geometry.

Resistive MHDFluid ionDKIGKI
0.0086ωA 0.0098ωA 0.002ωA 0.0022ωA 
Resistive MHDFluid ionDKIGKI
0.0086ωA 0.0098ωA 0.002ωA 0.0022ωA 

The growth rate of the tearing mode decreases for both drift kinetic and gyrokinetic ions. The finite Larmor radius (FLR) effects of thermal ions will not affect the growth rate of the tearing mode for these simulation parameters. When kinetic thermal ions are considered in the simulation, ion acoustic wave will be excited and carry parts of the total free energy away from the tearing mode. As a consequence, the drive of tearing mode is weakened.30 Moreover, the δϕ mode structures of both (1, 1) and (2, 1) modes become narrower than those of the fluid simulations as shown in Figs. 7 and 8. The reduction of the mode width and the growth rate are consistent, and probably due to the interaction between ions and tearing mode since the wave-particle interaction could become important at the edge of the mode envelop, i.e., r ≤ 0.29R0 and r ≥ 0.4R0, where vti ≥ γ/k. Nonetheless, the mechanisms for the reduction of the growth rate by the kinetic ion effects need to be studied further.

GTC can simulate fusion plasmas in general tokamak geometry. Implementing our resistive fluid model in the tokamak geometry, we can investigate the behavior of tearing mode in the realistic tokamak geometry. As a verification of the code capability, we use the simplest concentric tokamak for the fluid simulation of the (2, 1) resistive tearing mode. Using the same parameter as those of the cylindrical geometry and keeping the toroidal magnetic field curvature and gradient, we obtain the mode structure of the (2, 1) tearing mode in Fig. 9 from the tokamak simulations.

FIG. 9.

Comparison of radial mode structure of δA|| and δϕ from GTC fluid simulations in cylindrical and toroidal geometry. The position of mode rational surface q = 2 is at r/R0 = 0.343.

FIG. 9.

Comparison of radial mode structure of δA|| and δϕ from GTC fluid simulations in cylindrical and toroidal geometry. The position of mode rational surface q = 2 is at r/R0 = 0.343.

Close modal

The growth rate of the (2, 1) tearing mode is lower in the tokamak than that in the cylinder. For example, when νei = 0.001Ωp, the growth rate of (2, 1) tearing mode is 0.0028ωA in torus verse 0.00365ωA in cylinder. The reduced growth rate in the toroidal geometry is probably due to the favorable average curvature,31 Fig. 10 shows that mode structure of (2, 1) tearing mode in the toroidal geometry has a typical in-out asymmetry, which is the result of the in-out asymmetry of the toroidal magnetic field line.

FIG. 10.

Poloidal mode structures of δA|| (panel (a)) and δϕ (panel (b)) of (2, 1) tearing mode from GTC fluid simulation in tokamak geometry.

FIG. 10.

Poloidal mode structures of δA|| (panel (a)) and δϕ (panel (b)) of (2, 1) tearing mode from GTC fluid simulation in tokamak geometry.

Close modal

Gyrokinetic simulation of resistive tearing modes has been developed and verified in the GTC. GTC linear simulations in the fluid limit of the kink-tearing and resistive tearing modes in the cylindrical geometry agree well with the one-dimensional resistive magnetohydrodynamic eigenvalue and initial value codes. Ion kinetic effects are found to reduce the radial width and growth rate of the tearing modes. GTC simulations of the resistive tearing modes in the toroidal geometry find that the toroidicity reduces the growth rates. In the future work, gyrokinetic simulation of the resistive tearing modes will be extended to the nonlinear regime. By treating the tearing mode physics, neoclassical transport, and microturbulence on the same footing, the global gyrokinetic capability developed in this work will be applied in the future work to first-principles simulations of neoclassical tearing modes, which limit the performance of the burning plasmas such as ITER experiments.

The first author thanks helpful discussion with Dr. J. Bao, Z. X. Wang, and Professor Y. Q. Liu, and thanks H. S. Xie for the original eigenvalue code. This work was supported by U.S. DOE SciDAC GSEP Center and Grant No. DE-FG02-07ER54916, China National Magnetic Confinement Fusion Science Program (Grant Nos. 2013GB111000 and 2013GB112005), and National Natural Science Foundation of China (Grant No. 11205109). This research used resources of the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory (DOE Contract No. DE-AC05-00OR22725) and the National Energy Research Scientific Computing Center (DOE Contract No. DE-AC02-05CH11231).

1.
T.
Gosling
,
R. M.
Skoug
,
D. J.
McComas
, and
C. W.
Smith
,
J. Geophys. Res.
110
,
A01107
, doi: (
2005
).
2.
H. U.
Frey
,
T. D.
Phan
,
S. A.
Fuselier
, and
S. B.
Mende
,
Nature
426
,
533
537
(
2003
).
3.
M.
Øieroset
,
P. E.
Sandholt
,
W. F.
Denig
, and
S. W. H.
Cowley
,
J. Geophys. Res.
102
,
11349
11362
, doi: (
1997
).
4.
F. C.
Schüller
,
Plasma Phys. Controlled Fusion
37
,
A135
(
1995
).
5.
T. C.
Hender
,
J. C.
Wesley
,
J.
Bialek
,
A.
Bondeson
,
A. H.
Boozer
,
R. J.
Buttery
,
A.
Garofalo
,
T. P.
Goodman
,
R. S.
Granetz
,
Y.
Gribov
,
O.
Gruber
,
M.
Gryaznevich
,
G.
Giruzzi
,
S.
Günter
,
N.
Hayashi
,
P.
Helander
,
C. C.
Hegna
,
D. F.
Howell
,
D. A.
Humphreys
,
G. T. A.
Huysmans
,
A. W.
Hyatt
,
A.
Isayama
,
S. C.
Jardin
,
Y.
Kawano
,
A.
Kellman
,
C.
Kessel
,
H. R.
Koslowski
,
R. J.
La Haye
,
E.
Lazzaro
,
Y. Q.
Liu
,
V.
Lukash
,
J.
Manickam
,
S.
Medvedev
,
V.
Mertens
,
S. V.
Mirnov
,
Y.
Nakamura
,
G.
Navratil
,
M.
Okabayashi
,
T.
Ozeki
,
R.
Paccagnella
,
G.
Pautasso
,
F.
Porcelli
,
V. D.
Pustovitov
,
V.
Riccardo
,
M.
Sato
,
O.
Sauter
,
M. J.
Schaffer
,
M.
Shimada
,
P.
Sonato
,
E. J.
Strait
,
M.
Sugihara
,
M.
Takechi
,
A. D.
Turnbull
,
E.
Westerhof
,
D. G.
Whyte
,
R.
Yoshino
,
H.
Zohm
, and
the ITPA
MHD
, Disruption and Magnetic Control Topical Group.
Nucl. Fusion
47
,
6
(
2007
).
6.
H. P.
Furth
,
J.
Killeen
, and
M. N.
Rosenbluth
,
Phys. Fluids
6
,
459
(
1963
).
7.
H. P.
Furth
,
P. H.
Rutherford
, and
H.
Selberg
,
Phys. Fluids
16
,
1054
(
1973
).
8.
D.
Liu
and
L.
Chen
,
Plasma Phys. Controlled Fusion
53
,
062002
(
2011
).
9.
Z.-X.
Wang
,
L.
Wei
,
X.
Wang
,
S.
Zheng
, and
Y.
Liu
,
Nucl. Fusion
51
,
033003
(
2011
).
10.
H.
Cai
and
G.
Fu
,
Phys. Plasmas
19
,
072506
(
2012
).
11.
D. P.
Brennan
,
C. C.
Kim
, and
R. J.
La Haye
,
Nucl. Fusion
52
,
033004
(
2012
).
12.
H.
Cai
,
S.
Wang
,
Y.
Xu
,
J.
Cao
, and
D.
Li
,
Phys. Rev. Lett.
106
,
075002
(
2011
).
13.
Y.
Liu
,
R. J.
Hastie
, and
T. C.
Hender
,
Phys. Plasmas
19
,
092510
(
2012
).
14.
Z.
Lin
,
T. S.
Hahm
,
W. W.
Lee
,
W. M.
Tang
, and
R. B.
White
,
Science
281
,
1835
(
1998
).
15.
I.
Holod
,
W. L.
Zhang
,
Y.
Xiao
, and
Z.
Lin
,
Phys. Plasmas
16
,
122307
(
2009
).
16.
Z.
Lin
,
I.
Holod
,
L.
Chen
,
P. H.
Diamond
,
T. S.
Hahm
, and
S.
Ethier
,
Phys. Rev. Lett.
99
,
265003
(
2007
).
17.
Y.
Xiao
and
Z.
Lin
,
Phys. Rev. Lett.
103
,
085004
(
2009
).
18.
W.
Zhang
,
Z.
Lin
, and
L.
Chen
,
Phys. Rev. Lett.
101
,
095001
(
2008
).
19.
H. S.
Zhang
,
Z.
Lin
, and
I.
Holod
,
Phys. Rev. Lett.
109
,
025001
(
2012
).
20.
Z.
Wang
,
Z.
Lin
,
I.
Holod
,
W. W.
Heidbrink
, and
B.
Tobias
,
Phys. Rev. Lett.
111
,
145003
(
2013
).
21.
W.
Deng
,
Z.
Lin
, and
I.
Holod
,
Nucl. Fusion
52
,
023005
(
2012
).
22.
J.
McClenaghan
,
Z.
Lin
,
I.
Holod
, and
W.
Deng
,
Phys. Plasmas
21
,
122519
(
2014
).
23.
Z.
Lin
and
L.
Chen
,
Phys. Plasmas
8
,
1447
(
2001
).
24.
W. W.
Lee
,
J. Comput. Phys.
72
,
243
(
1987
).
25.
A. J.
Brizard
and
T. S.
Hahm
,
Rev. Mod. Phys.
79
,
421
(
2007
).
26.
R. B.
White
and
M. S.
Chance
,
Phys. Fluids
27
,
2455
(
1984
).
27.
A.
Mishchenko
and
A.
Zocco
,
Phys. Plasmas
19
,
122104
(
2012
).
28.
H.
Naitou
,
K.
Tsuda
,
W. W.
Lee
, and
R. D.
Sydora
,
Phys. Plasmas
2
,
4257
(
1995
).
29.
Y.
Nishimura
,
J. D.
Callen
, and
C. C.
Hegna
,
Phys. Plasmas
5
,
4292
(
1998
).
30.
M. N.
Bussac
,
D.
Edery
,
R.
Pellat
, and
J. L.
Soule
,
Phys. Rev. Lett.
40
,
1500
(
1978
).
31.
A. H.
Glasser
,
J. M.
Greene
, and
J. L.
Johnson
,
Phys. Fluids
18
,
875
(
1975
).