An analysis and preliminary experiment of the discharge characteristics of RF ion source with electromagnetic shielding

1.   Introduction

The hybrid kinetic–magnetohydrodynamic (MHD) simulation scheme [1], in comparison with the pure MHD and full particle-in-cell (PIC) simulation methods, has definite advantages that include the wave–particle energy exchange mechanism through the (inverse) Landau damping process and high efficiency in solving the evolution of background plasma. Especially, the hybrid kinetic–MHD simulation model has been widely adopted in simulations for studying the interaction between energetic particles (EPs) with various instabilities in tokamaks, such as the destabilizations of Alfvén eigenmodes and energetic particle modes (EPMs) by EPs through inverse Landau damping [2–4], and the influences of EPs on the linear stabilities of internal kink mode [5, 6] and resistive tearing modes (TMs) [7–11]. Over the past few decades, several hybrid codes have been developed, such as NIMROD [12], M3D-K [13], HMGC [14], MEGA [15], CLT-K [3], and the latest M3D-C1-K code [16].

It is well known that there are two coupling schemes in the hybrid kinetic–MHD model, i.e. pressure coupling and current coupling [1]. In the pressure coupling scheme, the perpendicular momentum evolution of EPs is ignored due to the assumption of n_h ≪ n (where n_h and n are the densities of EPs and background ions, respectively). In the current coupling scheme, all momentum components of EPs and background plasma are included self-consistently. In the previous version of the CLT-K code, the current coupling is adopted in the simulations [3, 4]. In this work, we bring the pressure coupling scheme into the CLT-K code. Numerical equivalences between these two coupling schemes have been strictly verified under different approximations in nonlinear simulations.

Using the hybrid simulation model, the influences of EPs on the linear stability of TM have been extensively studied through numerical simulations. However, contradictory results were reported in [8] and [11], especially for the situation of co-passing EPs. According to [8], the co-passing EPs play a stabilization role on the m/n ＝ 2/1 TM (where m and n represent the poloidal and toroidal mode numbers, respectively). But in [11], the co-passing EPs play a strong destabilization role on the m/n ＝ 2/1 TM. Based on numerical verifications of the equivalences of the current and pressure coupling schemes, the CLT-K code is adopted to reinvestigate the influences of the co-/counter-passing and trapped EPs on the linear stability of the m/n ＝ 2/1 TM. The co-passing and trapped EPs are found to stabilize the TM, but the counter-passing EPs have a destabilization effect on the TM. After exceeding the critical betas of EPs, the same branch of a high-frequency mode is excited by the co-/counter-passing and trapped EPs, which is identified as the m/n ＝ 2/1 EPM.

The outline of the present paper is as follows: section 2 introduces the simulation model used in the CLT-K code; section 3 presents the verification of the equivalences between the current and pressure coupling schemes; the simulation results about the influences of the co-/counter-passing and trapped EPs on the linear stability of the m/n ＝ 2/1 TM are presented in section 4; the resonant excitations of the m/n ＝ 2/1 EPM are discussed in section 5; finally, the summary and discussion of the simulation results are given in section 6.

2.   Simulation model

For background plasma, CLT-K solves the compressible MHD equations [3, 17–22], covering the plasma density $\rho ,$ thermal pressure $p,$ plasma velocity $\mathbf{v},$ magnetic field $\mathbf{B},$ electric field $\mathbf{E},$ and current density $\mathbf{J},$

with

All variables are nondimensionalized as $\mathbf{B}/B_{00}\to \mathbf{B},$ $\mathbf{x}/a\to \mathbf{x},$ $\rho /{\rho }_{00}\to \rho ,$ $\mathbf{v}/v_{\mathrm{A}}\to \mathbf{v},$ $t/{\tau }_{\mathrm{A}}\to t,$ $\gamma \mathrm{and}\omega /{\omega }_{\mathrm{A}}\to \gamma \mathrm{and}\omega ,$ $p/\left(B_{\mathrm{axis0}}^2/{\mu }_0\right)\to p,$ $\mathbf{J}/\left(B_{00}/{\mu }_{0}a\right)\to \mathbf{J},$ $\mathbf{E}/\left(v_{\mathrm{A}}{B}_{\mathrm{axis0}}\right)\to \mathbf{E},$ $\eta /\left({\mu }_0{v}_{\mathrm{A}}a\right)\to \eta ,$ and $D,\kappa ,\mathrm{and}\nu /\left(v_{\mathrm{A}}a\right)\to D,\kappa ,\mathrm{and}\nu ,$ where $a$ is the poloidal minor radius; $v_{\mathrm{A}}＝B_{\mathrm{axis0}}/\sqrt{{\mu }_0{\rho }_{\mathrm{axis0}}}$ is the Alfvén speed; and ${\tau }_{\mathrm{A}}＝a/v_{\mathrm{A}}$ and ${\omega }_{\mathrm{A}}＝v_{\mathrm{A}}/a$ are the Alfvén time and Alfvén frequency, respectively. The variable with subscript '$\mathrm{axis0}$' corresponds to its initial value at the magnetic axis. $\mathrm{\Gamma }$ is the adiabatic index (5/3).

For the MHD portion of the CLT-K code, i.e. the CLT code, the cylindrical coordinate system (R, φ, Z) is adopted, where R is the major radius, φ is the toroidal angle, and Z is along the vertical direction. The uniform and rectangular grid is used in the R–Z directions. For the discretization in space, the 4th-order finite difference method is employed in the R and Z directions, while in the φ direction either the 4th-order finite difference or pseudo-spectrum method can be used. As for the time advance, the 4th-order Runge–Kutta scheme is applied.

The contributions of EPs are taken into consideration in momentum equation with the EPs' current ${\mathbf{J}}_{\mathrm{h}}$ or pressure tensor ${\mathbf{P}}_{\mathrm{h}}$ in the hybrid kinetic–MHD model [1], i.e.

where the subscript 'h' denotes the EPs. The term $-Z_{\mathrm{h}}{n}_{\mathrm{h}}e{\mathbf{E}}_{\perp }$ (where $Z_{\mathrm{h}}$ is the charge number of EPs, $n_{\mathrm{h}}$ is the particle density of EPs, and $e$ is the unit charge) counteracts the $\mathbf{E}\times \mathbf{B}$ drift in the EP current and is eliminated in equation (7).

To push EPs, we adopt the guiding-center equations of motion [23], i.e.

where $M$ is the EP mass. The modified fields ${\mathbf{B}}^{☆}$ and ${\mathbf{E}}^{☆}$ are expressed in terms of the new effective electrostatic potential ${\mathrm{\Phi }}^{☆}$ and the vector potential ${\mathbf{A}}^{☆},$ i.e. ${\mathbf{B}}^{☆}＝\mathrm{▽}\times {\mathbf{A}}^{☆},$ $B_{||}^{☆}＝{\mathbf{B}}^{☆}·\mathbf{b},$ ($\mathbf{b}＝\mathbf{B}/B$), and ${\mathbf{E}}^{☆}＝-\mathrm{▽}{\mathrm{\Phi }}^{☆}-{\partial }_t{\mathbf{A}}^{☆}.$ ${\mathrm{\Phi }}^{☆}$ and ${\mathbf{A}}^{☆}$ are defined as

where $\mu ,$ $\mathrm{\Phi }$, and $\mathbf{A}$ are, respectively, the magnetic moment, electric potential, and magnetic vector potential. The parallel and perpendicular components are defined based on the unit vector of the magnetic field $\mathbf{b},$ e.g. ${\mathbf{v}}_{||}＝\left(\mathbf{v}·\mathbf{b}\right)\mathbf{b}$ and ${\mathbf{v}}_{\perp }＝\mathbf{b}\times \mathbf{v}\times \mathbf{b}.$

Ignoring polarization current [3, 4], ${\mathbf{J}}_{\mathrm{h}}$ in equation (7) mainly includes the current of guiding-center motion ${\mathbf{J}}_{\mathrm{GC}}$ and magnetization current ${\mathbf{J}}_{\mathrm{MAG}},$

where

The effective guiding-center drift velocities (curvature drift ${\mathbf{v}}_{\mathrm{curvature}},$ magnetic drift ${\mathbf{v}}_{\mathrm{▽}B},$ and parallel motion ${\mathbf{v}}_B$) are

By introducing the following intermediate variables, i.e. particle density $n_{\mathrm{h}},$ particle flux density $n_{\mathrm{h}}{V}_{\mathrm{h}||},$ perpendicular pressure component $P_{\mathrm{h}\perp }$, and parallel pressure component $P_{\mathrm{h}||},$

${\mathbf{J}}_{\mathrm{h}}$ in equation (13) can be expressed as

Meanwhile, ${\mathbf{P}}_{\mathrm{h}}$ is written in the Chew–Goldberger–Low (CGL) form [24]

The $\delta f$ method [3, 25, 26] is used in CLT-K to solve EP distribution function. At present, the contribution of EPs to equilibrium is also ignored.

3.   Equivalence of the two coupling schemes

In the hybrid kinetic–MHD model, the additional term introduced by EPs in the momentum equation is the ${\mathbf{J}}_{\mathrm{h}}\times \mathbf{B}$ in equation (7) or ${(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}})}_{\perp }$ in equation (8). According to equations (23) and (24), the analytical equivalence between the two additional terms can be easily proved as

The parallel components of ${\mathbf{J}}_{\mathrm{h}}$ and $\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}$ do not contribute to the momentum evolution of EPs. It can also be found that with the CGL pressure formulation as equation (24), the variation of the perpendicular momentum of the EPs is completely ignored due to the perfect cancellation between ${\mathbf{J}}_{\mathrm{h}}\times \mathbf{B}$ and ${(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}})}_{\perp }.$ Consequently, the current coupling and pressure coupling schemes for equations (7) and (8) are exactly the same when the full distribution function of EPs is used to calculate the total ${\mathbf{P}}_{\mathrm{h}}$ and ${\mathbf{J}}_{\mathrm{h}}$ as equations (23) and (24).

3.1   Perturbed forms of the momentum equations

The initial equilibria of the hybrid simulation are usually obtained by solving the Grad–Shafranov equation. However, the pressure obtained with this method without considering EPs is an isotropic scalar, which means the anisotropy of the EP pressure tensor is not included. The direct inclusion of EPs in simulation will result in an unbalance between the total Lorentz force and the pressures (tensor) of the background plasma and the EPs, i.e. $\mathcal{F}=\mathbf{J}_0 \times \mathbf{B}_0-\nabla p_0-\left(\nabla \cdot \mathbf{P}_{\mathrm{h} 0}\right)_{\perp \mathbf{b}_0} \neq 0$. (The subscript means perpendicular to the initial magnetic field.) This force unbalance will significantly affect the early simulation result, especially in linear stages. To avoid this issue, the initial residual force $\mathcal{F}$ is subtracted in equations (7) and (8).

For the current coupling scheme in equation (7), we subtract initial contributions of the background plasma, that is, $\mathcal{F}=\left(\mathbf{J}_0-\mathbf{J}_{\mathrm{h} 0}\right) \times \mathbf{B}_0-\nabla p_0 \neq 0$ Equation (7) becomes

where

For the pressure coupling scheme in equation (8), we subtract contributions of the total plasma, that is, $\mathcal{F}=\mathbf{J}_0 \times \mathbf{B}_0-\nabla p_0-\left(\nabla \cdot \mathbf{P}_{\mathrm{h} 0}\right)_{\perp \mathbf{b}_0} \neq 0$ Equation (8) becomes

where

and

If the initial force balance condition is satisfied in an equilibrium considering both contributions of the background plasmas and the EPs, i.e. ${\mathbf{J}}_0\times {\mathbf{B}}_0-\mathrm{▽}{p}_0-{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h0}}\right)}_{\perp {\mathbf{b}}_0}＝0$ and $\left({\mathbf{J}}_0-{\mathbf{J}}_{\mathrm{h0}}\right)\times {\mathbf{B}}_0-\mathrm{▽}{p}_0＝0,$ equations (26)–(30) will be completely equivalent to the total forms of the momentum equations of equations (7), (8), (23), and (24). Thus, there are two parts of contributions in the perturbed EP current and pressure tensor. The black perturbation terms without underline in equations (26)–(30) are caused by the variation of the EP distribution function $\delta f,$ while the underlined terms in equations (26)–(30) mainly result from the magnetic field perturbation $\delta \mathbf{B}.$ In the following part, it will be proved that the current and pressure coupling schemes are equivalent when only considering the $\delta f$ contribution or both terms resulted from $\delta f$ and $\delta \mathbf{B},$ respectively.

3.2   Equivalence with the reduced $\delta {\mathbf{P}}_{\mathrm{h}}$ and $\delta {\mathbf{J}}_{\mathrm{h}}:$ including only $\delta f$

First, we ignore the variation of magnetic field $\delta \mathbf{B},$ i.e. only the black terms without underline related to $\delta f$ in equations (26)–(30) are included. Then, equations (27) and (29) become

In addition, ${\mathbf{J}}_{\mathrm{h}0}\times \delta \mathbf{B}$ in equation (26) should also be ignored. The momentum equations become

and

It can be found that the equivalence of the current and pressure coupling schemes based on equations (31)–(34) can be proved similarly following equation (25). In comparison with equations (23) and (24), the only difference that appears in equations (31) and (32) is that the scalar variables of $n_{\mathrm{h}}{V}_{||},$ $P_{\mathrm{h}||},$ and $P_{\mathrm{h}\perp }$ are replaced by their variations as $\delta \left(n_{\mathrm{h}}{V}_{||}\right),$ $\delta P_{\mathrm{h}||},$ and $\delta P_{\mathrm{h}\perp }.$

For most hybrid kinetic–MHD codes, the simplified formulae of equations (31)–(34) are widely adopted. Representatively, equation (31) in the current coupled momentum equation is used in MEGA [15] and the early version of the CLT-K code [3], and equation (32) is used in the pressure tensor coupled momentum equation in NIMROD [12], M3D-K [13], and HMGC [14].

Recently, we have introduced the current and pressure coupling schemes in the CLT-K code [27]. With equations (31)–(34), we have conducted a set of fully nonlinear simulations for the n ＝ 1 TAE to demonstrate the numerical equivalence of these two coupling schemes. A mesh with 200 × 16 × 200 grids in ($R,\phi ,Z$) is adopted. About 40 markers are loaded in each grid cell. These parameters guarantee numerical convergence. The initial equilibrium and simulation parameters used are the same as those in previous work [27]. Here, we use a low resistivity as $\eta ＝{10}^{-7}.$ The nonlinear evolutions of TAE are shown in figure 1. The results obtained with the current and pressure coupling schemes are exactly the same for the n ＝ 1 TAE branch and the n ≠ 1 nonlinear sidebands. The eigenmode structure and frequency obtained from CLT-K can be found in figure 5 of [27] and will not be presented here again.

Figure  1.  Time evolutions of the kinetic energy in the n ＝ 1 TAE nonlinear simulations with the current coupling (solid lines) and pressure coupling (dotted lines) schemes. Only the variation of the EP distribution function $\delta f$ is considered.

DownLoad: Full-Size Img PowerPoint

3.3   Equivalence with the full $\delta {\mathbf{P}}_{\mathrm{h}}$ and $\delta {\mathbf{J}}_{\mathrm{h}}:$ including both $\delta f$ and $\delta \mathbf{B}$

Another more complicated situation is to include the contributions of both $\delta f$ and $\delta \mathbf{B}$ as equations (26)–(30). The equivalence of the current and pressure coupling schemes with equations (26)–(30) is hard to be demonstrated analytically. We adopted equations (26)–(30) in the latest CLT-K code and carried out a set of nonlinear simulations for the n ＝ 1 TAE with the same parameters, as shown in figure 1. According to the results shown in figure 2, the numerical equivalence between the two coupling schemes is well satisfied. In the linear stage of TAE, the results between the two methods, as shown in figures 1 and 2, are almost the same. However, the major difference made by the variation of the magnetic field appears significantly in the nonlinear saturation stage of TAE, especially for the n ＝ 3 and 4 branches. A possible explanation is that by including the underlined terms related to the magnetic field perturbation in $\delta {\mathbf{J}}_{\mathrm{h}}\times \mathbf{B}$ and $\delta \left[{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }\right],$ the nonlinearity of the magnetic field will be taken into account self-consistently. With the underlined terms in equations (26)–(30), the nonlinearity of the magnetic field is of the order of $O\left(\delta B^2\right).$ Otherwise, the nonlinearity included in the reduced $\delta {\mathbf{J}}_{\mathrm{h}}\times \mathbf{B}$ or $\delta \left[{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }\right]$ (without the underlined terms) mainly comes from the perturbed EP distribution function $\delta f$ and the magnetic field perturbation $\delta \mathbf{B},$ i.e. $O\left(\delta f\delta B\right).$ The nonlinearity of the magnetic field tends to generate nonlinear branches with higher n, e.g. n ＝ 3 and 4. On the other hand, in our simulations, we only considered the n ＝ 1 component of $\delta f.$ In consequence, the nonlinearity of $\delta f$ and $\delta \mathbf{B},$ i.e. $O\left(\delta f\delta B\right),$ will mainly generate the lower-n nonlinear branches, e.g. n ＝ 0, 1, and 2.

Figure  2.  Time evolutions of the kinetic energy in the n ＝ 1 TAE nonlinear simulations with the current coupling (solid lines) and pressure coupling (dotted lines) schemes. Both the variations of the EP distribution function $\delta f$ and the magnetic field $\delta \mathbf{B}$ are considered.

DownLoad: Full-Size Img PowerPoint

3.4   Brief discussion

Table 1 briefly summarizes the equations and the features of three forms of the current or pressure coupling schemes, i.e. (1) complete form; (2) perturbed form, including terms contributed from both $\delta f$ and $\delta \mathbf{B};$ and (3) perturbed form mainly considering the contribution of $\delta f$. For each form, the current and pressure coupling schemes are basically equivalent to each other.

Table  1.  Comparison of different forms of the current or pressure coupling schemes.

Items	Complete form	Perturbed form, including and	Perturbed form, including only
Current coupling	Equations (7) and (23)	Equations (26) and (27)	Equations (31) and (33)
Pressure coupling	Equations (8) and (24)	Equations (28)–(30)	Equations (32) and (34)
Additional simplifications	None	The initial force unbalance $\mathcal{F}$ caused by EPs is ignored	The underlined terms in equations (26)–(30) related to and the initial force unbalance $\mathcal{F}$ caused by EPs are ignored
TAE results	Not available	Larger saturation amplitudes of high-n branches	Smaller saturation amplitudes of high-n branches

 | Show Table

DownLoad: CSV

The complete forms of the coupling equations are strictly derived from the momentum equations of EPs and the background plasma. However, these complete coupling equations are unfriendly to the present equilibrium that cannot handle the anisotropic contribution of EPs. As a result, the initial force unbalance $\left[\mathcal{F}=\mathbf{J}_0 \times \mathbf{B}_0-\nabla p_0-\left(\nabla \cdot \mathbf{P}_{\mathrm{h} 0}\right)_{\perp \mathbf{b}_0} \neq 0\right]$ will cause a large perturbation at the initial stage of the simulation.

After we subtracted the initial force unbalance $\mathcal{F}$ in the momentum equations, we obtained the first perturbed form of the coupling equations. This form of equations is exactly the same as the complete form when the initial equilibrium condition, including the contribution of EPs, is satisfied, i.e. $\mathcal{F}=\mathbf{J}_0 \times \mathbf{B}_0-\nabla p_0-\left(\nabla \cdot \mathbf{P}_{\mathrm{h} 0}\right)_{\perp \mathbf{b}_0}=0$ The expressions of and are complicated because they include perturbations from both of EPs and the perturbed magnetic field

In general, we focus on wave–particle resonance problems, especially for the linear stage of modes. Then, the contribution from $\delta f$ will be more important. In this case, it is reasonable to ignore the underlined perturbed terms associated with $\delta \mathbf{B}$, and we can get the second simplified perturbed form of the coupling equations. The last perturbed form of the coupling equations is also widely adopted by various hybrid codes, as listed in the introduction and in section 3.2.

Simulation results from these two forms of perturbed coupling equations are rather similar for TAE, especially for the linear stage. The major difference is found in the saturation levels of high-n branches, which might be resulted from different nonlinear terms included in the equations, as discussed in section 3.3.

4.   Influences of EPs on the stability of the m/n ＝ 2/1 TM

In this section, we will systematically discuss the influences of EPs with different distribution functions on the linear stability of the m/n ＝ 2/1 TM. Furthermore, to understand the results, we will analyze the perturbed potential energy $\delta W_{\mathrm{h}}$ caused by EPs. At TM's linear stage, because the magnetic field perturbation $\delta \mathbf{B}$ is small, we adopt the simplified coupling schemes (only including the contribution of $\delta f$). The results from the current and pressure coupling schemes are exactly the same.

4.1   Simulation parameters and the initial equilibrium

We use a mesh with 200 × 16 × 200 grids in ($R,\phi ,Z$) for all cases, and about 40 markers are loaded in each grid cell. These parameters guarantee the numerical convergence. The anisotropic slowing down distribution with the pitch angle $\mathrm{\Lambda }＝\mu B_0/E$ is used for EPs

with 0 (passing EPs) and 1 (trapped EPs), and is the poloidal flux being averaged over the EP orbit [3], where and are, respectively, the poloidal fluxes at the magnetic axis and the boundary. The EP's Larmor radius $\varrho_{\mathrm{h}}$ with the perpendicular speed of is 0.07 Nevertheless, the finite Larmor radius effect is not included in CLT-K because we employ the drift-kinetic model to describe the guiding-center orbits of EPs. The of EPs is scanned with its value at the magnetic axis varying from 0.154% to 3% while the EP pressure profile is fixed. The uniform resistivity used in the simulations is fixed to be Other dissipation coefficients are

Here, we adopt a zero beta (${\beta }_{\mathrm{b}}$ ＝ 0) tokamak equilibrium, which is obtained from the QSOLVER code [28] with shifted circular flux surfaces. The aspect ratio is $R_0/a$ ＝ 5.76. The plasma density is uniform. The profile of the safety factor is plotted in figure 3(a). The m/n ＝ 2/1 TM's instability parameter $\mathrm{\Delta }\text{'}$ is 14.8 [29], which means the m/n ＝ 2/1 TM is unstable. The saturated mode structure $E_{\phi }$ of the m/n ＝ 2/1 TM without EPs is plotted in figure 3(b). We mainly discuss the influences of EPs on the m/n ＝ 2/1 TM. The n ＝ 1 perturbed $\delta {\mathbf{J}}_{\mathrm{h}}$ and $\delta {\mathbf{P}}_{\mathrm{h}}$ are included with only the contribution of $\delta f$ as equations (31)–(34) in the simulations. In addition, $\delta f$ associated with the EP response can be separated into the adiabatic and non-adiabatic parts [8, 30]. The adiabatic term $-\xi ·\mathrm{▽}{f}_0$ in $\delta f$ is the fluidlike response of the EPs (where $\xi$ is the plasma displacement), which will be subtracted in $\delta f$ to investigate its role in the (de)stabilization of TMs. The non-adiabatic part of $\delta f$ as $\delta f+\xi ·\mathrm{▽}{f}_0,$ mainly accounts for the resonance and the energy exchange between EPs and waves. To get comparable results, the initial EP distribution and tokamak equilibrium are the same as those adopted by Cai et al [8].

Figure  3.  (a) Initial profile of equilibrium's safety factor $q.$ (b) The mode structure of the m/n ＝ 2/1 TM without EPs in the nonlinear saturated stage. The initial q ＝ 2 surface is denoted by the dashed line.

DownLoad: Full-Size Img PowerPoint

4.2   Influences of EPs on the linear stability of TM

4.2.1   Co-passing EPs

We first present the simulation results of the m/n ＝ 2/1 TM in the presence of co-passing EPs. The linear growth rate of the m/n ＝ 2/1 TM versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of co-passing EPs is plotted in figure 4(a). When including the total response of EPs, as shown by red circles in figure 4(a), the growth rate of the TM decreases slightly with the increase of ${\beta }_{\mathrm{h}}^{\mathrm{c}},$ which means that EPs, in general, play a weak stabilization effect on the TM. The frequency of the TM is very close to the analytical diamagnetic drift frequency ${\omega }_{\mathrm{h}}^{*}$ of EPs at the q ＝ 2 surface; see figure 4(b). The negative frequency indicates that the mode rotates in the direction of the ion-diamagnetic drift. The mode frequency in figure 4(b) is missing when ${\beta }_{\mathrm{h}}^{\mathrm{c}}>1.0\%.$ This is because a high-frequency m/n ＝ 2/1 EPM is excited, which results in difficulty in the identification of the TM frequency (if TM exists).

Figure  4.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 mode versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of co-passing EPs. The results with (without) the adiabatic term are indicated by red circles (blue squares). The diamagnetic drift frequency of the analytical EPs at the q ＝ 2 surface is also plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

Next, we subtract the adiabatic contribution in $\delta f$ of EPs, that is, only the non-adiabatic response of EPs is included in simulations. As shown by blue squares in figure 4(a), without the adiabatic term, the non-adiabatic response of co-passing EPs plays a strong effective role in the stabilization of the TM. For 0.40% $\le {\beta }_{\mathrm{h}}^{\mathrm{c}}\le$ 1.0%, the TM is fully stabilized by the non-adiabatic response of EPs. The mode frequency can still be measured because we initially add an m/n ＝ 2/1 TM-like perturbation around the q ＝ 2 rational surface to trigger the fast growth of the m/n ＝ 2/1 TM. With the adiabatic response of co-passing EPs, the growth rate of TM is much larger than that without the adiabatic response. Consequently, the adiabatic response of co-passing EPs has a strong destabilization effect on the TM. For ${\beta }_{\mathrm{h}}^{\mathrm{c}}<0.40\%,$ the mode frequencies of the TM without adiabatic term are also quite close to those in the presence of both the adiabatic and non-adiabatic terms.

4.2.2   Counter-passing EPs

The influence of counter-passing EPs on the linear stability of the m/n ＝ 2/1 TM is also investigated. As shown in figure 5, within the scanned range of ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ ($<1.15\%$), the dominant mode in the system is the TM with low diamagnetic drift frequency. The total response of counter-passing EPs has a small destabilization effect on the TM, as indicated by red circles in figure 5(a). Meanwhile, the non-adiabatic response of counter-passing EPs significantly destabilizes the TM with the growth rate even increased by a factor of 3. Thus, the adiabatic response of counter-passing EPs greatly stabilizes the TM. The frequency of the TM resulting from counter-passing EPs is also comparable with the diamagnetic drift frequency of the analytical EPs, as shown in figure 5(b).

Figure  5.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 TM versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of counter-passing EPs. The results with (without) the adiabatic term are indicated by red circles (blue squares). The diamagnetic drift frequency of the analytical EPs at the q ＝ 2 surface is also plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

4.2.3   Trapped EPs

Here, we briefly study the role of trapped EPs () on the linear stability of the m/n ＝ 2/1 TM. Cai et al adopt a large value of () [8], which results in a quasi-isotropic EP distribution instead of a trapped one. However, here, we adopt a narrow distribution function in pitch angel space () so that most of EPs belong to trapped particles. The results from CLT-K are plotted in figure 6. Qualitatively, the trapped EPs have a similar influence on the stability of the TM as the co-passing EPs. Overall, the total response of trapped EPs slightly stabilizes the TM, but the non-adiabatic response dramatically stabilizes the TM. Nevertheless, the stabilization effect from the non-adiabatic response of trapped EPs is strong so that the TM is completely stabilized when $\beta_{\mathrm{h}}^{\mathrm{c}} \gtrsim 0.3 \%$.

Figure  6.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 TM versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of trapped EPs. The results with (without) the adiabatic term are indicated by red circles (blue squares). The diamagnetic drift frequency of the analytical EPs at the q ＝ 2 surface is also plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

Similar to the situations of passing EPs, trapped EPs also result in a rotation of the TM in the frequency of , as shown in figure 6(b). However, when including only the non-adiabatic response of trapped EPs and $\beta_{\mathrm{h}}^{\mathrm{c}} \gtrsim 0.8 \%$, the mode frequency disappears. This is because the mode structure of the TM is unobservable, which indicates the full stabilization of the TM by trapped EPs.

4.3   Mechanisms of EPs on the linear stability of the TM

The linear results about passing EPs on the stability of the m/n ＝ 2/1 TM by CLT-K support the analytical and numerical results of Cai et al [7, 8]. For trapped EPs, the distribution function in our simulations greatly differs from that of [8], so the comparison is of lesser significance. According to the theoretical work by Cai et al [7], the non-adiabatic responses of passing EPs mainly influence the perturbed parallel current $\delta J_{||}$ in the ideal region, which results in the modifications of $\mathrm{\Delta }\text{'}$ and the linear stability of the TM. In our hybrid simulation, the numerical noise caused by PIC simulation makes $\delta J_{||}$ difficult to be analyzed. In the following part, to understand the influence of EPs on the stability of the TM, we will analyze the potential energy $\delta W_{\mathrm{h}}$ caused by EPs with different distribution functions. By taking the pressure coupling as an example, $\delta W_{\mathrm{h}}$ can be written as

where $\mathit{\kappa }＝\mathbf{b}·\mathrm{▽}\mathbf{b}$ is the magnetic field curvature. The same formula of $\delta W_{\mathrm{h}}$ can be obtained by replacing ${\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }$ in equation (36) by ${\mathbf{J}}_{\mathrm{h}}\times \mathbf{B}$ for the current coupling scheme.

We first analyze the adiabatic response of EPs. For simulation cases in the presence of total response of EPs, the poloidal distribution of the integral term $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}^{\mathrm{adiabatic}}\right)}_{\perp }$ at the $\phi ＝0$ cross-section is plotted in figure 7. Co-passing and trapped EPs have stronger adiabatic responses in the low field side (LFS), while the main response of counter-passing EPs appears in the high field side (HFS). The adiabatic response of EPs peaks at the X and O points of the TM, where the plasma displacement $\xi$ is the strongest. The poloidal asymmetry of $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}^{\mathrm{adiabatic}}\right)}_{\perp }$ is mainly caused by the orbit shifts of EPs. The (de)stabilization roles of the adiabatic responses of the EPs can be understood based on the sign of $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}^{\mathrm{adiabatic}}\right)}_{\perp }.$ Specifically, in figure 7(a), the adiabatic response of co-passing EPs is negative at the X and O points of the TM, which indicates a destabilization effect on the TM. On the contrary, the adiabatic response of counter-passing EPs mainly stabilizes the TM since $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}^{\mathrm{adiabatic}}\right)}_{\perp }$ is positive in figure 7(b). For trapped EPs in figure 7(c), the adiabatic response dominantly plays a destabilization role in the LFS but a weak stabilization role in the HFS. Overall, the adiabatic response of trapped EPs destabilizes the TM. The adiabatic responses of EPs are well consistent with the linear simulation results shown in figures 4–6.

Figure  7.  The 2D distribution of $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}^{\mathrm{adiabatic}}\right)}_{\perp }$ in the presence of total response of EPs (${\beta }_{\mathrm{h}}^{\mathrm{c}}＝0.57\%;$ $t＝2991.2{\tau }_{\mathrm{A}}$): (a) co-passing EPs, (b) counter-passing EPs, and (c) trapped EPs. The dashed line indicates the q ＝ 2 surface. The Poincaré plots of magnetic field lines are presented by dotted lines. Note that the perturbed magnetic fields are enlarged to clearly show the locations of the X and O points of magnetic islands.

DownLoad: Full-Size Img PowerPoint

Next, we analyze $\delta W_{\mathrm{h}}$ generated by the total response of EPs so as to understand the behind (de)stabilization mechanisms of the TM by EPs. The integral term $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }$ contributed by the total $\delta f$ contains high numerical noise. So, we integrate $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }$ along the toroidal direction and plot the 2D distribution of ${\displaystyle \frac{1}{2}}{\displaystyle \int \xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }R\mathrm{d}\phi }$ in figure 8. The value of $\delta W_{\mathrm{h}}$ is also calculated by performing 3D volume integral of $\xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }.$ The total response of co- and counter-passing EPs dominates in the LFS and HFS, respectively, while the response of trapped EPs has a quasi-symmetric poloidal distribution. The different distributions of $\delta W_{\mathrm{h}}$ are related to the various features of EP orbits. The orbits of co- (counter)-passing EPs are significantly shifted toward LFS (HFS). The m/n ＝ 2/1 TM is located around the q ＝ 2 rational surface. As a result, the dominant interaction of co- (counter)-passing EPs and TM will naturally take place in the LFS (HFS) of the q ＝ 2 rational surface. However, for the case with trapped EPs, we used a moderate pitch-angle parameter ${\mathrm{\Lambda }}_0＝1.0,$ and particles with a large pitch angle will also be randomly loaded in the HFS initially, which together result in a non-negligible portion of barely passing EPs in addition to the trapped EPs. In our simulations, it is difficult to separate or remove the barely passing EPs since their pitch angle is also large. For trapped EPs, they are mainly located in the LFS due to their orbits, which dominantly contributes to $\delta W_{\mathrm{h}}$ in the LFS. But, for the barely passing EPs, they will spend much longer time due to their small $v_{||}$ in the HFS. As a result, the interaction of the barely passing EPs and the m/n ＝ 2/1 TM will be obvious in the HFS. Thus, the distribution of $\delta W_{\mathrm{h}}$ for the case with trapped EPs is quasi-symmetric in the poloidal section. On the other hand, inside (outside) the q ＝ 2 surface, EPs with different distributions play stabilizing (destabilizing) roles on the TM. The antisymmetric distribution of ${\displaystyle \frac{1}{2}}{\displaystyle \int \xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }R\mathrm{d}\phi }$ around the q ＝ 2 surface is related to the odd symmetric mode structure of the m/n ＝ 2/1 TM. The integrated values of $\delta W_{\mathrm{h}}$ in figure 8 show that co-passing and trapped EPs play stabilizing roles on the TM ($\delta W_{\mathrm{h}}>0$), while counter-passing EPs destabilize the TM ($\delta W_{\mathrm{h}}<0$). The analyses of $\delta W_{\mathrm{h}}$ also agree with the linear results of figures 4–6.

Figure  8.  The 2D distribution of ${\displaystyle \frac{1}{2}}{\displaystyle \int \xi ·{\left(\mathrm{▽}·{\mathbf{P}}_{\mathrm{h}}\right)}_{\perp }R\mathrm{d}\phi }$ in the presence of the total response of EPs (${\beta }_{\mathrm{h}}^{\mathrm{c}}＝0.57\%;$ $t＝2991.2{\tau }_{\mathrm{A}}$): (a) co-passing EPs, (b) counter-passing EPs, and (c) trapped EPs. The integrated values of $\delta W_{\mathrm{h}}$ are marked out. The dashed line indicates the q ＝ 2 surface.

DownLoad: Full-Size Img PowerPoint

5.   Resonant excitation of the m/n ＝ 2/1 EPM

In this section, we will briefly discuss the excitation of the m/n ＝ 2/1 EPM by EPs with different distribution functions. The simulation parameters are the same as in section 4, but the ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of EPs are increased to higher values.

For co-passing EPs, the growth rates and frequencies of the m/n ＝ 2/1 mode are plotted in figure 9. When the dominant mode is the low-frequency $\left(\omega \simeq \omega_{\mathrm{h}}^*\right) m / n=2 / 1$ TM. However, when exceeds the threshold an m/n ＝ 2/1 mode with $\omega \simeq-0.025 \omega_{\mathrm{A}}$ is excited by co-passing EPs. The frequency of this m/n ＝ 2/1 mode is much larger than and clearly intersects with the m ＝ 2 component of the shear Alfvén wave (SAW) continuum, as shown in figure 10(a). Correspondingly, the mode structure of the m/n ＝ 2/1 poloidal electric field is plotted in figure 10(b). The mode structure of this high-frequency mode extends toward the magnetic axis and is totally different from that of the m/n ＝ 2/1 TM, as shown in figure 3(b). Consequently, this high-frequency mode is identified as the m/n ＝ 2/1 EPM. In figure 9(a), with the increasing the growth rate of the TM first decreases and then that of EPM increases obviously.

Figure  9.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 mode versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of co-passing EPs. The results with (without) the adiabatic term are indicated with red circles (blue squares). The diamagnetic drift frequency of the analytical EPs is plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

Figure  10.  (a) Mode-frequency spectrum and (b) mode structure of the m/n ＝ 2/1 poloidal electric field $E_{\theta }$ in the presence of co-passing EPs with ${\beta }_{\mathrm{h}}^{\mathrm{c}}＝1.43\%.$ The n ＝ 1 SAW continuum is plotted with black solid lines (a). The dashed line in (b) indicates the q ＝ 2 surface.

DownLoad: Full-Size Img PowerPoint

The same branch of EPM (with similar mode frequencies and mode structures as in figure 10) is also observed in cases with counter-passing and trapped EPs. Respectively shown in figures 11 and 12, the excitation thresholds of the m/n ＝ 2/1 EPM are about $2\%$ for counter-passing EPs and $1.5\%$ for trapped EPs. The frequency spectra of the EPM in the presence of counter-passing and trapped EPs are very similar to those in figure 10. Meanwhile, as shown in figures 9 and 11, the adiabatic responses of co- and counter-passing EPs have very limited influences on the excitation of EPM. However, the adiabatic response of trapped EPs greatly promotes the excitation of EPM, which is similar to the result reported in [8]. This should be caused by the wide banana orbit of trapped EPs so that the trapped EPs have a much larger adiabatic response on EPM than passing EPs. In the recent theoretical work by Zhang et al [31], the excitation of the m/n ＝ 2/1 EPM by the non-adiabatic response of trapped EPs has been systematically discussed. Nonetheless, for passing EPs, they assumed the orbit width is small, so that the leading order of non-adiabatic response from passing EPs can be ignored and the m/n ＝ 2/1 EPM is hard to be excited [31]. According to our simulations, co-passing EPs have the smallest threshold for the excitation of the EPM, indicating that the finite orbit width effect of passing EPs is non-negligible.

Figure  11.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 mode versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of counter-passing EPs. The results with (without) the adiabatic term are indicated by red circles (blue squares). The diamagnetic drift frequency of the analytical EPs is plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

Figure  12.  (a) Linear growth rates and (b) mode frequencies (at the q ＝ 2 surface) of the m/n ＝ 2/1 mode versus ${\beta }_{\mathrm{h}}^{\mathrm{c}}$ of trapped EPs. The results with (without) the adiabatic term are indicated by red circles (blue squares). The diamagnetic drift frequency of the analytical EPs is plotted by green stars in (b).

DownLoad: Full-Size Img PowerPoint

In the end, we present resonant conditions of EPs with the m/n ＝ 2/1 EPM. The resonance condition can be written as $\omega +p{\omega }_{\theta }-n{\omega }_{\phi }＝0,$ where ${\omega }_{\theta }$ and ${\omega }_{\phi }$ are the poloidal and toroidal orbit frequencies, respectively; $p$ is an integer with $p＝l+m,$ while $l$ is the Fourier component of the poloidal motion of EPs. The $\delta f$ structures of EPs with different distribution functions are plotted in figure 13. For co-passing EPs ($\mathrm{\Lambda }＝0.25$) in figure 13(a), the dominant resonance condition is $p＝1.$ The $p＝1$ resonance condition is different from previous report [32], where the $p＝2$ resonance condition was found by M3D-K simulations. This difference should be caused by the different energies of co-passing EPs adopted in simulations of CLT-K ($v_0＝0.5v_{\mathrm{A}}$) and M3D-K ($v_0>v_{\mathrm{A}}$). On the other hand, with our simulation parameters and initial equilibrium, p ＝ 2, n ＝ 1, $\omega ＝-0.020{\omega }_{\mathrm{A}}$ resonance condition of co-passing EPs, and the m/n ＝ 2/1 EPM is not covered by the phase space of EPs. For counter-passing EPs, the $p＝3$ resonance condition is found to be dominant, which also results in a significant radial transport of counter-passing EPs; see figure 13(b). The large increase of $\delta f$ near the plasma boundary, i.e. the upper boundary in figure 13(b), indicates a significant loss of counter-passing EPs caused by the EPM. Meanwhile, the precession frequency of trapped EPs matches the EPM frequency. Thus, the dominant resonance condition is $p＝0$ for trapped EPs. Also, higher resonance conditions are observed in figure 13(c), such as $p＝1$ and 2.

Figure  13.  The $\delta f$ structures of EPs for cases in the presence of (a) co-passing ($\mathrm{\Lambda }＝0.25$), (b) counter-passing ($\mathrm{\Lambda }＝0.25$), and (c) trapped ($\mathrm{\Lambda }＝1.0$) EPs. The different resonance conditions between EPs and TM are plotted with different lines.

DownLoad: Full-Size Img PowerPoint

6.   Conclusion and discussion

In this paper, we strictly verified the equivalences of the current and pressure coupling schemes under different approximations. For most of hybrid MHD–kinetic simulations, only the contribution of $\delta f$ is considered, which captures the most important wave–particle resonance and the adiabatic response of EPs. However, if the magnetic perturbation caused by instabilities grows to a large amplitude, the contribution of $\delta \mathbf{B}$ in the perturbed pressure tensor or the current of EPs should also be retained. The significant difference by the $\delta \mathbf{B}$ contribution can be found in the nonlinear stage of TAE, as discussed in section 3. To describe the evolution of plasma in the presence of EPs more accurately, a self-consistent equilibrium is necessary to include the initial contribution of EPs properly, e.g. the anisotropic pressure equilibrium, including parallel and perpendicular pressure components. The optimization of such an equilibrium is complicated but worth pursuing further, because it directly influences the widely used Grad–Shafranov model and the closure of the MHD equations.

Then, influences of co-/counter-passing and trapped EPs on the linear stability of the m/n ＝ 2/1 TM are carefully reinvestigated by adopting the first approximation method in CLT-K, i.e. considering only the $\delta f$ contribution of EPs. The contribution of $\delta \mathbf{B}$ in the perturbed pressure tensor or the current of EPs is ignored because of the small amplitude of $\delta \mathbf{B}$ in the linear stage. The CLT-K simulations show that non-adiabatic responses of co-passing and trapped EPs stabilize the TM, while that of counter-passing EPs has the opposite effect. Moreover, the adiabatic responses of EPs with different distributions tend to counteract the non-adiabatic (de)stabilization effects of co-/counter-passing and trapped EPs. The behind (de)stabilization mechanisms are numerically analyzed for the first time based on the perturbed potential energy $\delta W_{\mathrm{h}}$ produced by EPs. The simulation results and numerical analyses both present a consistent conclusion as Cai et al [7, 8] for co-/counter-passing EPs.

Finally, we studied the general excitations of the high-frequency m/n ＝ 2/1 EPM by EPs. The transition of TM toward EPM occurs when the beta of EPs is large enough to overcome damping mechanisms, e.g. continuum damping. In our simulations, the beta thresholds for exciting the EPM by different distributed EPs are in the order as follows: the co-passing EPs have the smallest beta threshold, then the trapped EPs in second, and the beta threshold of counter-passing EPs is the largest. The simulation results partially confirm the analytical research by Zhang et al [31], which focused on the trapped EPs. In addition, the excitation behaviors of EPM by co-/counter-passing EPs via different resonance conditions [32] are also systematically discussed.