Verification of gyrokinetic particle simulation of current-driven instability in fusion plasmas. II. Resistive tearing mode

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 sun¹ 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.⁴ 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 ITER⁵ 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 D_e 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.⁶ Some of the simulation work of the tearing mode have been performed on large parallel simulation codes, such as M3D¹⁰ and NIMROD.¹¹ 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,¹⁸ and Alfvén eigenmodes.^19,20 The implementation of the equilibrium current²¹ enables simulation of current-driven instabilities in the toroidal geometry with kinetic effects, such as the kink mode²² 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,²³ by assuming ω/k_‖v_te ≪ 1; here, ω is the mode frequency of interest, k_‖ is the parallel wave number, and v_te 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:

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(fe−fe0)$ 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

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

where

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=∇×δA∥b0$ 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 electron²¹ 

and the parallel momentum equation

Here, the total electron density $ne$⁠, parallel flow $u∥e$⁠, and pressure $pe$ are the sum of their equilibrium and perturbed parts, i.e., $ne=ne0+δne$⁠, $u∥e=u∥e0+δu∥e$⁠, 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

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

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

The electrostatic potential is calculated from the gyrokinetic Poisson's equation²⁴ 

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

In these equations, the ion guiding center density $ni$ and current $u∥i$ can be calculated from the standard gyrokinetic model.²⁵ 

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 field²⁶ can be represented as either covariant form

or contravariant form

and the Jacobian is

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.²¹ 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 $ρi2∇⊥2≪1$ becomes

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

Together with the massless electron parallel momentum equation

and parallel Ampère's law

where $∇⊥2=1/r(∂(r∂/∂r)−m2/r2)$⁠, $∇||=(n−m/q)BT/RB0$⁠, and $δB=∇×δA∥b0$⁠. 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/R₀ = 0.5, major radius R₀ = 100 cm, magnetic field B₀ = 10 000G, equilibrium electron density on magnetic axis n_e0 = 10¹⁴/cm³, and equilibrium electron beta β_e = $pe0/(8πB02)$ = 0.004. Note that the position of q = 1 mode rational surface is r = 0.25R₀. For ideal MHD, in which $∂δA||/∂t=−cb0⋅∇δϕ$⁠, we get the ideal kink mode structures in Fig. 1.

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/R₀ = 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.1R₀, which is much larger than the thermal ion gyroradius 0.0022R₀ 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.

Furth et al.⁶ 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.6155−0.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.

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 structure^7,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.

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.

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 T_i = T_e = 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.

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.³⁰ 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.29R₀ and r ≥ 0.4R₀, where v_ti ≥ γ/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.

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,³¹ 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.

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).