Global simulations of kinetic-magnetohydrodynamic processes with energetic electrons in tokamak plasmas

Considering the smallness of realistic electron mass, we apply the drift-kinetic model for EE species that induces the necessary Landau resonance while ignores the FLR effect. The linearized drift-kinetic equation in five dimensional phase space uses guiding center position R, magnetic moment $\mu = m_ev_\perp^2/2B_0$ and parallel velocity $v_{||}$ as independent variables, which reads

$$\begin{align} \begin{aligned} L_0\delta f_{h} + \delta L^L f_{h0} = 0 \end{aligned}, \end{align} \tag{ 1 }$$

where $\delta f_{h}$ and $f_{h0}$ are the EE perturbed and equilibrium distributions, L₀ and $\delta L^L$ are the equilibrium and linear perturbed propagators given by

$$\begin{align*} \begin{aligned} L_0 = \frac{\partial }{\partial t} + \left(v_{||}\mathbf{b_0} + \mathbf{v_d}\right)\cdot\nabla - \frac{\mu}{m_e B_0}\mathbf{B_0^*}\cdot\nabla B_0 \frac{\partial }{\partial v_{||}} \end{aligned} \end{align*}$$

and

$$\begin{align*} \begin{aligned} \delta L^L & = \left(v_{||}\frac{\mathbf{\delta B}}{B_0}+\mathbf{v_E}\right)\cdot\nabla - \left[\frac{\mu}{m_eB_0}\mathbf{\delta B}\cdot\nabla B_0 \right.\\ & \quad \left. +\ \frac{q_e}{m_e}\left(\frac{\mathbf{B_0^*}}{B_0}\cdot\nabla\delta \phi + \frac{1}{c}\frac{\partial\delta A_{||}}{\partial t}\right)\right]\frac{\partial }{\partial v_{||}} \end{aligned}, \end{align*}$$

q_e and m_e are electron charge and mass, δφ and $\delta A_{||}$ are the electrostatic potential and parallel vector potential, $\mathbf{b_0} = \mathbf{B_0}/B_0$ is the unit vector of the equilibrium magnetic field, $\mathbf{B_0^*} = \mathbf{B_0} + \left(B_0v_{||}/\Omega_{ce}\right)\nabla\times\mathbf{b_0}$, $\Omega_{ce} = q_eB_0/cm_e$ is the electron cyclotron frequency, $\mathbf{\delta B} = \nabla\times \left(\delta A_{||}\mathbf{b_0}\right)$ represents the perturbed magnetic field, and $\mathbf{v_d}$ and $\mathbf{v_E}$ are the magnetic drift and $\mathbf{E}\times \mathbf{B}$ drift, which read

$$\begin{align*} \begin{aligned} \mathbf{v_d} = \frac{v_{||}^2}{\Omega_{ce}}\nabla\times\mathbf{b_0} + \frac{\mu\mathbf{b_0}\times\nabla B_0}{m_e\Omega_{ce}} \end{aligned} \end{align*}$$

and

$$\begin{align*} \begin{aligned} \mathbf{v_E} = \frac{c\mathbf{b_0}\times\nabla\delta\phi}{B_0} \end{aligned}. \end{align*}$$

For EE-driven Alfvenic instabilities, the linear unstable spectra can be influenced by the EE velocity distribution relying on the specific auxiliary heating methods such as ECRH and LHCD [13]. Though the EE velocity distribution is not unique in experiments, we define effective EE temperatures for analysis convenience, namely, $T_{||h0} = \int m_ev_{||}^2f_{h0} \mathbf{dv}/n_{h0}$ and $T_{\perp h0} = \int \mu B_0f_{h0} \mathbf{dv}/n_{h0}$, where $n_{h0} = \int f_{h0} \mathbf{dv}$ is the EE density. In current work that illustrates the model scheme, we utilize the isotropic Maxwellian to approximate the EE equilibrium distribution, i.e. $f_{h0} = n_{h0} (\frac{m_e}{2\pi T_{h0}})^{3/2}exp (-\frac{m_ev_{||}^2+2\mu B_0}{2T_{h0}})$ with $T_{h0} = T_{||h0} = T_{\perp h0}$, and other EE velocity distributions such as slowing-down will be investigated in future work. We further separate the perturbed distribution $\delta f_h$ in equation (1) into the adiabatic part and non-adiabatic part

$$\begin{align} \begin{aligned} \delta f_{h} = \delta f^A + \delta K \end{aligned}. \end{align} \tag{ 2 }$$

The adiabatic distribution is determined by the terms related to fast parallel dynamics in equation (1), which are in the leading order of $\omega/(k_{||}v_{||})$ as

$$\begin{align} \begin{aligned} & v_{||}\mathbf{b_0}\cdot\nabla\delta f^A + v_{||}\frac{1}{B_0}\mathbf{\delta B}\cdot\nabla f_{h0}\Big{\lvert}_{v_\perp}\\ & \quad - \frac{q_e}{m_e}\left(\mathbf{b_0}\cdot\nabla\delta\phi + \frac{1}{c}\frac{\partial \delta A_{||}}{\partial t}\right)\frac{\partial f_{h0}}{\partial v_{||}} = 0 \end{aligned}, \end{align} \tag{ 3 }$$

where $\nabla f_{h0}{\lvert}_{v_\perp} = (\nabla n_{h0}/n_{h0})f_{h0} + [(m_ev_{||}^2 + 2\mu B_0)/(2T_{h0})-3/2](\nabla T_{h0}/T_{h0})f_{h0}$ and $\partial f_{h0}/\partial v_{||} = -(m_ev_{||}/T_{h0})f_{h0}$ for Maxwellian $f_{h0}$. We define a new variable δψ according to

$$\begin{align} \begin{aligned} \frac{\partial \delta A_{||}}{\partial t} = -c\mathbf{b_0}\cdot\nabla\delta \psi \end{aligned}, \end{align} \tag{ 4 }$$

and $\delta\psi = \omega\delta A_{||}/(ck_{||})$ is adopted for following derivations with convenience. Substituting equation (4) into equation (3) and using the ansatz $\partial_t = -i\omega$ and $\mathbf{b_0}\cdot\nabla = ik_{||}$, $\delta f^A$ can be readily solved from equation (3)

$$\begin{align} \begin{aligned} \delta f^A & = - \frac{q_e}{T_{h0}}\left(\delta\phi - \delta \psi\right)f_{h0} \\ & \quad -\frac{q_e}{T_{h0}}\delta \psi\left[\frac{\omega_{*n,h}}{\omega}+\left(\frac{m_ev_{||}^2 + 2\mu B_0}{2T_{h0}}-\frac{3}{2}\right)\frac{\omega_{*T,h}}{\omega}\right]f_{h0} \end{aligned}, \end{align} \tag{ 5 }$$

where $\omega_{*n,h} = -i\frac{cT_{h0}}{q_eB_0}\mathbf{b_0}\times\frac{\nabla n_{h0}}{n_{h0}}\cdot\nabla$ and $\omega_{*T,h} = -i\frac{c}{q_eB_0}\mathbf{b_0}\times\nabla T_{h0}\cdot\nabla$. Note that the form of $\delta f^A$ in equation (5) is also adopted in gyrokinetic theory [27, 28] and simulation [29], which contains both the adiabatic response to the parallel electric field and convective response to the perturbed magnetic field.

From equations (1), (2) and (5), the governing equation for the non-adiabatic distribution $\delta K$ can be written as

$$\begin{align} \begin{aligned} &\left[-i\left(\omega - \omega_{d}\right) + v_{||}\mathbf{b_0}\cdot\nabla - \frac{\mu}{m_e B_0}\mathbf{B_0^*}\cdot\nabla B_0 \frac{\partial }{\partial v_{||}}\right]\delta K\\ & \quad = - i\frac{q_e}{T_{h0}}\omega \left(1-\frac{\omega_{*p,h}}{\omega}\right)\left(\delta\phi - \delta\psi\right)f_{h0}\\ & \qquad -i\frac{q_e}{T_{h0}}\omega_{d}\left(1-\frac{\omega_{*p,h}}{\omega}\right)\delta\psi f_{h0} \end{aligned}, \end{align} \tag{ 6 }$$

where $\omega_{d} = -i\mathbf{v_d}\cdot\nabla$ is the magnetic drift frequency and $\omega_{*p,h} = \omega_{*n,h} + [(m_ev_{||}^2 + 2\mu B_0)/(2T_{h0})-3/2]\omega_{*T,h}$ is the diamagnetic frequency. It is also noted that $C_s/V_A = \sqrt{\beta_e/2}$, $v_\mathrm{the}/V_A = \sqrt{\beta_em_i/(2m_e)}$ and $C_s/v_\mathrm{the} = \sqrt{m_e/m_i}$, where C_s , V_A and $v_\mathrm{the}$ denote ion sound speed, Alfven speed and electron thermal speed, respectively, $\beta_e = 8\pi n_eT_e/B_0^2$ and $m_e/m_i\ll\beta_e\ll 1$ in tokamak plasmas. With these orderings, one can solve $\delta K$ from equation (6) for passing and trapped EE particles, respectively.

First, in the regime of $k_{||}v_{||}\gg\omega\sim\omega_{d}$ for passing EE particles, equation (6) reduces to

$$\begin{align} \begin{aligned} v_{||}\mathbf{b_0}\cdot\nabla\delta K^{p} & = - i\frac{q_e}{T_{h0}}\omega \left(1-\frac{\omega_{*p,h}}{\omega}\right)\left(\delta\phi - \delta\psi\right)f_{h0}\\ & \quad -i\frac{q_e}{T_{h0}}\omega_{d}\left(1-\frac{\omega_{*p,h}}{\omega}\right)\delta\psi f_{h0} \end{aligned}, \end{align} \tag{ 7 }$$

where we have $\delta K^{p}\sim 0$ since $v_{||}\mathbf{b_0}\cdot\nabla\to\infty$, and $\delta K^{p}$ contribution to perturbed density and pressure can be ignored. However, $v_{||}\delta K^{p}$ is finite and gives rise to a fraction of parallel current density carried by passing EEs.

Second, EE particles interact and exchange energy with MHD modes primarily through the precessional drift resonance, while the bounce and transit frequencies are too high to resonate with typical AEs, as demonstrated and confirmed by recent theories and simulations [14, 16, 17, 19]. According to these characteristics, we have $\omega\sim\overline{\omega}_{d}\ll\omega_{b}$ for trapped EE particles, where $\omega_{b} = 2\pi/\tau_{b}$ and $\tau_{b} = \oint \left(\mathrm{d}l/v_{||}\right)$ is the bounce period, and $\overline{\omega}_{d} = \oint \omega_{d}\left(\mathrm{d}l/v_{||}\right)/\tau_{b}$ is the precession frequency, and l is the traveled distance of trapped particle along magnetic field line in one bounce motion period. Defining $\omega_{b}\partial/\partial\eta = v_{||}\mathbf{b_0}\cdot\nabla$ for trapped particle with η being the bounce angle coordinate and $\mathrm{d}\eta = \omega_{b}\left(\mathrm{d}l/v_{||}\right)$, we can transform equation (6) into banana orbit center frame using $\delta K^t = \delta K^t_b \mathrm{exp}\left(i\alpha\right)$ with $\alpha = -\int \mathrm{d}\eta\left(\omega_{d}-\overline{\omega}_{d}\right)/\omega_{b}$ [30], and then perform bounce average $\overline{\left(\cdots\right)} = \oint \left(\cdots\right)\left(\mathrm{d}l/v_{||}\right)/\tau_{b}$, i.e. $\overline{\left(\cdots\right)} = \oint \left(\cdots\right)\mathrm{d}\eta/2\pi$. Note that $\omega_{d}-\overline{\omega}_{d}\ll\omega_{b}$, we have $\mathrm{exp}\left(i\alpha\right)\approx1$ which means that $\delta K^t \simeq \delta K^t_b$ and finite orbit width (FOW) can be dropped for EE particles as a higher order effect, then the solution of trapped EE non-adiabatic response is

$$\begin{align} \begin{aligned} \delta K^t \simeq \delta K^t_b & = \underbrace{\frac{\omega}{\omega-\overline{\omega}_{d}}\frac{q_e}{T_{h0}}\left(1-\frac{\omega_{*p,h}}{\omega}\right)\left(\overline{\delta\phi} -\overline{\delta\psi}\right)f_{h0}}_{\left\{I\right\}}\\ & \quad \underbrace{+ \frac{1}{\omega-\overline{\omega}_{d}}\frac{q_e}{T_{h0}}\left(1-\frac{\omega_{*p,h}}{\omega}\right)\overline{\omega_{d}\delta\psi} f_{h0} }_{\left\{II\right\}} \end{aligned}, \end{align} \tag{ 8 }$$

where term {I} corresponds to the non-ideal MHD effect $\Delta \phi = \delta\phi-\delta\psi$, and term {II} drives the bad-curvature instabilities with a minimum threshold $\omega_{*p,h}\gt\omega$. For most instabilities in tokamak characterized with 'flute-like' mode structures, δφ and δψ peak around the rational surfaces where $|nq - m|\ll 1$, then we can further simplify equation (8) using $\overline{\delta\phi} \approx \delta\phi$, $\overline{\delta\psi} \approx \delta\psi$, and $\overline{\omega_{d}\delta\psi}\approx\overline{\omega}_{d}\delta\psi$. The more general validity regime of these simplifications is $\mathrm{exp}\left[-i\left(m-nq\right)\theta\right]\approx 1$ [16], which only requires $|\theta|\ll 1/|nq-m|$. Thus, for instabilities with moderate $|nq-m|\unicode{x2A7D}1$ that deviate from rational surface, one can reduce trapped EE fraction by decreasing the θ domain of banana orbit mirror throat to more deeply-trapped regime, which still hold $\overline{\delta\phi} \approx \delta\phi$, $\overline{\delta\psi} \approx \delta\psi$, and $\overline{\omega_{d}\delta\psi}\approx\overline{\omega}_{d}\delta\psi$ for model accuracy. (Note that these simplifications are applied for deriving the EE moments in section 2.2).

To couple drift-kinetic EEs with the existing Landau-fluid model of bulk plasmas, one needs to derive the EE continuity equation by integrating equation (1) in velocity space using $\langle\cdots\rangle_v = \int \mathbf{dv} = \frac{2\pi B_0}{m_e}\int \mathrm{d}v_{||}\mathrm{d}\mu$, i.e.

$$\begin{align} \begin{aligned} \frac{\partial \delta n_{h}}{\partial t} &+ \mathbf{v_E}\cdot\nabla n_{h0} + n_{h0}\mathbf{B_0}\cdot\nabla\left(\frac{\delta u_{||h}}{B_0}\right) + 2cn_{h0}\nabla\delta\phi \cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0}\\ & + \frac{c}{q_e} \nabla\left(\delta P_{||h}^A + \delta P_{\perp h}^A\right)\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0} + \left\langle \boldsymbol{v_d}\cdot\nabla\delta K\right\rangle_v = 0, \end{aligned} \end{align} \tag{ 9 }$$

where

$$\begin{align} \delta n_{h} & = \int \left(\delta f^A + \delta K\right) \mathbf{dv}, \end{align} \tag{ 10 }$$

$$\begin{align} n_{h0}\delta u_{||h} & = \int v_{||}\left(\delta f^A + \delta K\right) \mathbf{dv}, \end{align} \tag{ 11 }$$

$$\begin{align} \delta P_{||h}^A & = \int m_ev_{||}^2\delta f^A \mathbf{dv} = - q_en_{h0}\left(\delta\phi - \delta \psi\right)\nonumber\\ & \quad - q_en_{h0} \left(\frac{\omega_{*n,h}}{\omega}+\frac{\omega_{*T,h}}{\omega}\right)\delta \psi \end{align} \tag{ 12 }$$

and

$$\begin{align} \begin{aligned} \delta P_{\perp h}^A = \int \mu B_0\delta f^A \mathbf{dv} & = - q_en_{h0}\left(\delta\phi - \delta \psi\right)\\ &\quad -q_en_{h0}\left(\frac{\omega_{*n,h}}{\omega}+\frac{\omega_{*T,h}}{\omega}\right)\delta \psi \end{aligned}. \end{align} \tag{ 13 }$$

Taking the moment of equation (5), the adiabatic components in equations (10) and (11) are obtained as

$$\begin{align} \begin{aligned} \delta n_{h}^A = \int \delta f^A \mathbf{dv} = -\frac{q_en_{h0}}{T_{h0}}\left(\delta\phi - \delta \psi\right) - \frac{q_en_{h0}}{T_{h0}}\frac{\omega_{*n,h}}{\omega}\delta\psi \end{aligned} \end{align} \tag{ 14 }$$

and

$$\begin{align} \begin{aligned} \delta u_{||h}^A = \frac{1}{n_{h0}}\int v_{||}\delta f^A \mathbf{dv} = 0 \end{aligned}, \end{align} \tag{ 15 }$$

Equations (12)–(15) are referred to as adiabatic EE moments, where the superscript 'A' stands for adiabatic.

Then we shall derive the non-adiabatic EE moments. Since the passing EEs move much faster than Alfven speed of which response is almost adiabatic to shear Alfven wave and ion sound wave, $\delta K^p$ in equation (7) does not contribute to density and pressure perturbations but can lead to a finite current in high-$v_{||}$ regime. However, the precessional drift of trapped EE can be close to or below the Alfven speed, and $\delta K^t$ in equation (8) are responsible for the non-adiabatic EE density and pressure. Moreover, it has been shown that the deeply-trapped EEs can effectively interact with most drift-type instabilities through precessional drift resonance [14, 17, 19], thus we apply the deeply-trapped approximation for computing the precession frequency in equation (8)

$$\begin{align} \begin{aligned} \overline{\omega}_d \approx \omega_{d0} \end{aligned}, \end{align} \tag{ 16 }$$

where $\omega_{d0} = \omega_d\Big\lvert_{\theta = 0}$ is the magnetic drift frequency at the outer midplane which only contains the normal curvature related term, while the geodesic curvature related term vanishes in the bounce-averaged dynamics and does not contribute to the precessional drift resonance [31]. The form of $\omega_{d0}$ in general flux coordinates is given by equation (69), and the validation of deeply-trapped approximation on precession frequency is also shown in section 4. Based on the deeply-trapped approximation, the last term in equation (9) can be written as

$$\begin{align} \begin{aligned} \left\langle \boldsymbol{v_d}\cdot\nabla\delta K\right\rangle_v \simeq\left\langle \boldsymbol{v_d}\cdot\nabla\delta K^t\right\rangle_v = i\frac{2}{T_{h0}}\omega_{D0,h}\delta P_{h}^{NA} \end{aligned} \end{align} \tag{ 17 }$$

where $\omega_{D0,h} = - i \frac{cT_{h0}}{q_eB_0}\mathbf{b_0}\times\boldsymbol{\kappa}\cdot\nabla\Big\lvert_{\theta = 0}$, and $\delta P_{h}^{NA} \simeq \frac{1}{2}\int_\mathrm{trap} \delta K^t E\mathbf{dv}$ represents the non-adiabatic EE pressure. We use pitch angle $\lambda = \mu B_a/\left(m_eE\right)$ and energy per unit mass $E = v^2/2$ as the two dimensions $\left(\lambda,E\right)$ of velocity space instead of $\left(v_{||}, \mu\right)$, where B_a is the on-axis magnetic field strength, then the velocity space integration with trapped particles becomes

$$\begin{align} \begin{aligned} \int_\mathrm{trap} \mathbf{dv} = \sqrt{2}\pi\frac{B_0}{B_a}\sum_{\sigma}\int_{\lambda_\mathrm{low}}^{B_a/B_0}\frac{1}{\sqrt{\left(1-\lambda\frac{B_0}{B_a}\right)}}\mathrm{d}\lambda\int_{0}^{+\infty}\sqrt{E}\mathrm{d}E \end{aligned}, \end{align} \tag{ 18 }$$

where $\sigma = \mathrm{sign}(v_{||})$, and $B_a/B_\mathrm{max}\unicode{x2A7D}\lambda_\mathrm{low} \unicode{x2A7D} B_a/B_0$ is the lower cutoff of trapped particle pitch angle. For Maxwellian equilibrium distribution $f_{h0} = n_{h0}\left(\frac{m_e}{2\pi T_{h0}}\right)^{3/2}exp\left(-\frac{m_eE}{T_{h0}}\right)$, the trapped particle fraction f_t at B₀ location can be expressed as

$$\begin{align} \begin{aligned} f_t = \frac{1}{n_{h0}}\int_\mathrm{trap} f_{h0}\mathbf{dv} = \sqrt{1-\lambda_\mathrm{low}\frac{B_0}{B_a}} \end{aligned}. \end{align} \tag{ 19 }$$

By counting all trapped particles at B₀ location with $\lambda_\mathrm{low} = B_a/B_\mathrm{max}$, it is seen that $f_t = 0$ in the high field side with $B_0 = B_\mathrm{max}$ and $f_t = \sqrt{1-B_\mathrm{min}/B_\mathrm{max}}$ in the low field side with $B_0 = B_\mathrm{min}$. Compared to former gyrokinetic theory that applies $f_t \approx \sqrt{2\epsilon}$ (where $\epsilon = r/R_0$) and ignores poloidal variation [31], equation (19) has both radial and poloidal dependences and reflects the realistic trapped particle distribution in the poloidal plane. In general, equation (18) can be simplified when the integrands do not rely on λ

$$\begin{align} \begin{aligned} \int_\mathrm{trap}\mathbf{dv} = 4\sqrt{2}\pi f_t\int_{0}^{+\infty}\sqrt{E}\mathrm{d}E \end{aligned}. \end{align} \tag{ 20 }$$

Thus, with the assumption in equation (16), it is straightforward to integrate $\delta K^t$ in equation (8) in velocity space using equation (20), and the non-adiabatic density and pressure perturbations are derived explicitly

$$\begin{align} \begin{aligned} \delta n_{h}^{NA} & = \int\delta K\mathbf{dv}\\ & = -f_t\frac{q_en_{h0}}{T_{h0}} \left[2\left(1-\frac{\omega_{*n,h}}{\omega}+\frac{3}{2}\frac{\omega_{*T,h}}{\omega}\right)\zeta R_1\left(\sqrt{\zeta}\right)\right.\\ & \quad \left. -\ 2\frac{\omega_{*T,h}}{\omega}\zeta R_3\left(\sqrt{\zeta}\right)\right]\left(\delta\phi - \delta\psi\right)\\ & \quad -f_t \frac{q_en_{h0}}{T_{h0}} \left[2\left(1-\frac{\omega_{*n,h}}{\omega}+\frac{3}{2}\frac{\omega_{*T,h}}{\omega}\right)R_3\left(\sqrt{\zeta}\right)\right.\\ & \quad \left.-\ 2\frac{\omega_{*T,h}}{\omega}R_5\left(\sqrt{\zeta}\right)\right]\delta\psi \end{aligned} \end{align} \tag{ 21 }$$

and

$$\begin{align} \begin{aligned} \delta P_{h}^{NA} & = \frac{1}{2}\int \delta K E\mathbf{dv}\\ & = -f_tq_en_{h0} \left[\left(1-\frac{\omega_{*n,h}}{\omega}+\frac{3}{2}\frac{\omega_{*T,h}}{\omega}\right)\zeta R_3\left(\sqrt{\zeta}\right)\right.\\ & \quad \left.-\ \frac{\omega_{*T,h}}{\omega}\zeta R_5\left(\sqrt{\zeta}\right)\right]\left(\delta\phi - \delta\psi\right) \\ & \quad -f_tq_en_{h0} \left[\left(1-\frac{\omega_{*n,h}}{\omega}+\frac{3}{2}\frac{\omega_{*T,h}}{\omega}\right)R_5\left(\sqrt{\zeta}\right)\right.\\ & \quad \left. -\ \frac{\omega_{*T,h}}{\omega}R_7\left(\sqrt{\zeta}\right)\right]\delta\psi \end{aligned}, \end{align} \tag{ 22 }$$

where $\zeta = \omega/\omega_{D0,h}$, the response functions are given by

$$\begin{align*} \begin{aligned} &R_1 \left(\sqrt{\zeta}\right) = 1+\sqrt{\zeta}Z\left(\sqrt{\zeta}\right),\\ &R_3 \left(\sqrt{\zeta}\right) = \frac{1}{2}+\zeta+\left(\zeta\right)^{3/2}Z\left(\sqrt{\zeta}\right),\\ &R_5 \left(\sqrt{\zeta}\right) = \frac{3}{4}+\frac{1}{2}\zeta+\zeta^2+\left(\zeta\right)^{5/2}Z\left(\sqrt{\zeta}\right),\\ &R_7\left(\sqrt{\zeta}\right) = \frac{15}{8}+\frac{3}{4}\zeta+\frac{1}{2}\zeta^2+\zeta^3+\left(\zeta\right)^{7/2}Z\left(\sqrt{\zeta}\right), \end{aligned} \end{align*}$$

and $Z\left(x\right) = \frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}\frac{\mathrm{exp}\left(-t^2\right)}{t-x}dt$ is the plasma dispersion function. Substituting equations (12)–(15), (17), (21) and (22) into equation (9), the non-adiabatic parallel velocity $\delta u_{||h}^{NA}$ is readily solved as

$$\begin{align} \begin{aligned} \delta u_{||h}^{NA} & = -\frac{q_e}{T_{h0}}\frac{\omega}{k_{||}}\left(1-f_t\right)\left(1-\frac{\omega_{*n,h}}{\omega}\right)\left(\delta\phi - \delta\psi\right) -2\frac{q_e}{T_{h0}}\frac{\omega}{k_{||}}\\ & \quad \times \left(\frac{\omega_{D,h}}{\omega}-f_t\frac{3}{4}\frac{\omega_{D0,h}}{\omega}\right)\left(1-\frac{\omega_{*n,h}}{\omega}-\frac{\omega_{*T,h}}{\omega}\right)\delta\psi \end{aligned}, \end{align} \tag{ 23 }$$

where $\omega_{D,h} = - i \frac{cT_{h0}}{q_eB_0}\mathbf{b_0}\times\boldsymbol{\kappa}\cdot\nabla$. In the derivation of equation (23), the resonance terms associated with $Z\left(\sqrt{\zeta}\right)$ in equations (21) and (22) exactly cancel out with each other through equation (9), which again proves that the trapped EE particles do not carry parallel current as we assumed before. Note that $\delta u_{||h}^{NA}$ can also be computed by integrating equation (7) in velocity space with passing EE fraction, which is consistent with equation (23) in the leading order.

Next, we couple the EE moments described by equations (12)–(15) and (21)–(23) to the bulk plasma Landau-fluid model in [25], and then yield the final fluid-kinetic hybrid model as follows. It should be pointed out that we use the electron charge symbol 'q_e ' in this work, while [25] uses the elementary charge symbol 'e', which lead to the sign differences of some terms since $q_e = -e$.

$$\begin{align} &\left[\frac{\partial}{\partial t}\left(1+ 0.75\rho_i^2\nabla_\perp^2\right)+ i\omega_{*p,i}\right]\frac{c}{V_A^2}\nabla_\perp^2\delta\phi\nonumber\\ & \quad +\mathbf{B_0}\cdot\nabla\left(\frac{1}{B_0}\nabla_\perp^2\delta A_{||}\right) - \frac{4\pi}{c}\boldsymbol{\delta B}\cdot\nabla\left(\frac{J_{||0}}{B_0}\right)\nonumber\\ & \quad -8\pi\left(\nabla\delta P_i+\nabla\delta P_e \right)\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0} \underbrace{-8\pi\nabla\delta P_{h}^{A}\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0}}_\mathrm{\{EE-IC\}}\nonumber\\ & \quad \underbrace{- i \frac{8\pi q_e}{cT_{h0}}\omega_{D0,h}\delta P_{h}^{NA}}_\mathrm{\{EE-KPC\}} = 0, \end{align} \tag{ 24 }$$

$$\begin{align} & \frac{\partial\delta A_{||}}{\partial t}\! =\! -c\mathbf{b_0}\cdot\nabla\delta\phi \!-\! \frac{cT_{e0}}{q_en_{e0}}\mathbf{b_0}\cdot\nabla\delta n_e \!-\! \frac{cT_{e0}}{q_en_{e0}B_0}\boldsymbol{\delta B}\cdot\nabla n_{e0}\nonumber\\ & \qquad\qquad - \frac{cm_e}{q_e}\sqrt{\frac{\pi}{2}}v_{the}|k_{||}|\delta u_{||e} + \frac{c^2}{4\pi}\eta_{||}\nabla_\perp^2\delta A_{||}, \end{align} \tag{ 25 }$$

$$\begin{align} & \frac{\partial \delta P_i}{\partial t} + \frac{c\mathbf{b_0}\times\nabla\delta\phi}{B_0}\cdot\nabla P_{i0} +2\Gamma_{i\perp} P_{i0}c\nabla\delta\phi\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0}\nonumber\\ & \quad + \Gamma_{i||} P_{i0}\mathbf{B_0}\cdot\nabla\left(\frac{\delta u_{||i}}{B_0}\right) -i\Gamma_{i\perp} \omega_{*p,i}Z_in_{i0}\rho_i^2\nabla_\perp^2\delta \phi\nonumber\\ & \quad + 2\Gamma_{i\perp} P_{i0}\frac{c}{Z_i}\nabla\delta T_i\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0} +2\Gamma_{i\perp} T_{i0}\frac{c}{Z_i}\nabla\delta P_i\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0}\nonumber\\ & \quad + n_{i0}\frac{2}{\sqrt{\pi}}\sqrt{2}v_{thi}|k_{||}|\delta T_i = 0, \end{align} \tag{ 26 }$$

$$\begin{align} & m_in_{i0}\frac{\partial\delta u_{||i}}{\partial t} = \frac{Z_in_{i0}}{q_en_{e0}}\left(\mathbf{b_0}\cdot\nabla\delta P_e + \frac{1}{B_0}\boldsymbol{\delta B}\cdot\nabla P_{e0}\right)\nonumber\\ & \qquad\qquad \quad - \left(\mathbf{b_0}\cdot\nabla\delta P_i + \frac{1}{B_0} \boldsymbol{\delta B}\cdot\nabla P_{i0}\right)\\ & \qquad\qquad \quad + Z_in_{i0}\frac{m_e}{q_e}\sqrt{\frac{\pi}{2}}v_{the}|k_{||}|\delta u_{||e} - Z_in_{i0}\frac{c}{4\pi}\eta_{||}\nabla_\perp^2\delta A_{||}, \end{align} \tag{ 27 }$$

$$\begin{align} & \frac{\partial \delta n_i}{\partial t} + \frac{c\mathbf{b_0}\times\nabla\delta\phi}{B_0}\cdot\nabla n_{i0} + 2cn_{i0}\nabla\delta\phi \cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0}\nonumber\\ & \quad + n_{i0}\mathbf{B_0}\cdot\nabla\left(\frac{\delta u_{||i}}{B_0}\right) - i \omega_{*p,i}\frac{Z_in_{i0}}{T_{i0}}\rho_i^2\nabla_\perp^2\delta\phi \\ & \quad + \frac{2c}{Z_i}\nabla\delta P_i\cdot\frac{\mathbf{b_0}\times\boldsymbol{\kappa}}{B_0} = 0, \end{align} \tag{ 28 }$$

where $\delta P_h^A = \delta P_{||h}^A = \delta P_{\perp h}^A$ in equation (24), $Z_in_{i0} + q_e\left(n_{e0}+n_{h0}\right) = 0$, and the definitions of bulk plasma variables are consistent with [25]. The two-moment and three-moment Landau closures are applied to thermal electrons in equation (25) and thermal ions in equation (26) respectively, which show good agreements with drift-kinetic theory on the response functions in the regime of $k_{||}v_{thi}\ll\omega\ll k_{||}v_{the}$ [32]. To close above equation set, we also have the equations for $\delta P_e$, $\delta T_e$ and $\delta T_i$ as

$$\begin{align} \begin{aligned} \delta P_e = \delta n_eT_{e0} + n_{e0}\delta T_e \end{aligned}, \end{align} \tag{ 29 }$$

$$\begin{align} \begin{aligned} \mathbf{b_0}\cdot\nabla\delta T_e + \frac{1}{B_0}\boldsymbol{\delta B}\cdot\nabla T_{e0} = 0 \end{aligned} \end{align} \tag{ 30 }$$

and

$$\begin{align} \begin{aligned} \delta T_i = \frac{1}{n_{i0}}\left(\delta P_i - \delta n_i T_{i0}\right) \end{aligned}. \end{align} \tag{ 31 }$$

Meanwhile, the thermal electron perturbed density $\delta n_e$ and parallel velocity $\delta u_{||e}$ in equations (24), (25) and (27) are calculated through the quasi-neutrality condition and parallel Ampere's law as

$$\begin{align} \begin{aligned} \delta n_e = -\frac{Z_i}{q_e}\delta n_i - \frac{c^2}{4\pi q_eV_A^2}\nabla_\perp^2\delta\phi \underbrace{- \left( \delta n_{h}^A + \delta n_{h}^{NA} \right)}_\mathrm{\{EE-density\}} \end{aligned} \end{align} \tag{ 32 }$$

and

$$\begin{align} \begin{aligned} q_en_{e0}\delta u_{||e} = -Z_in_{i0}\delta u_{||i} - \frac{c}{4\pi}\nabla_\perp^2\delta A_{||} \underbrace{-q_en_{h0}\delta u_{||h}}_\mathrm{\{EE-current\}} \end{aligned}. \end{align} \tag{ 33 }$$

Equations (4), (12)–(15), (21)–(23) and (24)–(33) form a closed system for the fluid-kinetic hybrid simulation model as shown in figure 1. The main characteristics are briefly summarized here: first the well-circulating and deeply-trapped approximations used for the integrals of EE moments do not break the EE continuity equation; second the EE and bulk plasma are coupled through the quasi-neutrality condition, parallel Ampere's law and vorticity equation in a non-perturbative manner (Note that non-perturbative here refers to solving the exact solution of nonlinear eigenmode equation in terms of ω, rather than full-f simulations of plasma nonlinear systems (i.e. nonlinearities in terms of distribution and field perturbations) without separating equilibrium and perturbed distributions); third in equation (24) the EE-KPC term (i.e. EE kinetic particle compression) is responsible for the dissipative excitation of AEs, and the EE-IC term (i.e. EE interchange from adiabatic fluid convection) contributes to the reactive MHD interchange drive.

Zoom In Zoom Out Reset image size

As claimed by equation (16) in section 2, we utilize the deeply-trapped approximation for calculating the precession frequency of trapped EE. To delineate its regime of validity, we compare equation (16) with the exact precession frequency by performing bounce average along the realistic banana orbit. The Boozer coordinates $\left(\psi,\theta,\zeta\right)$ is applied in MAS to describe the magnetic field (where $\psi$ is the poloidal magnetic flux, $\theta$ and $\zeta$ are the poloidal and toroidal Boozer angles respectively), which has the contravariant form $\mathbf{B_0} = q\left(\psi\right)\nabla\psi\times\nabla\theta - \nabla\psi\times\nabla\zeta$ and the covariant form $\mathbf{B_0} = I\left(\psi\right)\nabla\theta + g\left(\psi\right)\nabla\zeta$. Then the guiding-center equation of motion in equilibrium B-field can be expressed in $\left(\psi,\theta,\zeta\right)$ coordinates [33]

$$\begin{align} \dot{\psi} & = -\frac{c}{Z_\alpha}\left(\frac{m_\alpha v_{||}^2}{B_0} + \mu\right)\frac{1}{JB_0^2}\left(g\frac{\partial B_0}{\partial\theta}\right), \end{align} \tag{ 61 }$$

$$\begin{align} \dot{\theta} & = \frac{v_{||}}{JB_0} - \frac{cm_\alpha v_{||}^2}{Z_\alpha}\frac{1}{JB_0^2}\frac{\partial g}{\partial\psi} + \frac{c}{Z_\alpha}\left(\frac{m_\alpha v_{||}^2}{B_0}+\mu\right)\frac{1}{JB_0^2}g\frac{\partial B_0}{\partial\psi}, \end{align} \tag{ 62 }$$

$$\begin{align} \dot{\zeta} & = \frac{qv_{||}}{JB_0} + \frac{cm_\alpha v_{||}^2}{Z_\alpha}\frac{1}{JB_0^2}\frac{\partial I}{\partial\psi} -\frac{c}{Z_\alpha}\left(\frac{m_\alpha v_{||}^2}{B_0}+\mu\right)\frac{1}{JB_0^2}I\frac{\partial B_0}{\partial\psi} \end{align} \tag{ 63 }$$

and

$$\begin{align} \dot{\rho}_{||} = -\left(1-\rho_{||}g^{^{\prime}}\right)\frac{c}{Z_\alpha}\left(\frac{m_\alpha v_{||}^2}{B_0} + \mu\right)\frac{1}{JB_0^2}\frac{\partial B_0}{\partial\theta}, \end{align} \tag{ 64 }$$

where α denotes the particle species, $J = \left(\nabla\psi\times\nabla\theta\cdot\nabla\zeta\right)^{-1}$ is the Jacobian and $\rho_{||} = v_{||}/\Omega_{c\alpha}$. Given energy $E = 0.5v^2$ and pitch angle $\lambda = \mu B_a/\left(m_\alpha E\right)$ of a particle, parallel velocity $v_{||}$ can be written as

$$\begin{align} v_{||} = \pm\sqrt{2E\left(1-\lambda\frac{B_0}{B_a}\right)}, \end{align} \tag{ 65 }$$

and the bounce period can be calculated using equations (62) and (65) [33]

$$\begin{align} \tau_b = \oint \mathrm{d}\theta/\dot{\theta} \simeq \oint\frac{JB_0}{v_{||}}\mathrm{d}\theta. \end{align} \tag{ 66 }$$

Using equations (61)–(64) and (66), the exact precession frequency is then derived as [33]

$$\begin{align} \begin{aligned} \overline{\omega}_d & = n\overline{\mathbf{v_d}\cdot\nabla\left(\zeta-q\theta\right)}\\ & = \frac{n}{\tau_b}\left[-\oint \frac{1}{\dot{\theta}}\frac{c}{Z_\alpha}\left(\mu + \frac{m_\alpha v_{||}^2}{B_0}\right)\frac{\partial B_0}{\partial\psi}\mathrm{d}\theta \right.\\ & \left. \quad + \oint\frac{cm_\alpha v_{||}}{Z_\alpha B_0}\left(\frac{\partial I}{\partial\psi} + q\frac{\partial g}{\partial\psi}\right)\mathrm{d}\theta + \oint\rho_{||}g\frac{\mathrm{d}q}{\mathrm{d}\psi}\mathrm{d}\theta\right] \end{aligned}, \end{align} \tag{ 67 }$$

where n is the toroidal mode number, and the first term on the RHS is the leading order term that depends on both the radial ψ and poloidal θ coordinates, i.e. $\overline{\omega}_d = \overline{\omega}_d\left(\psi,\theta\right)$.

The magnetic drift frequency ω_d in Boozer coordinate can be expressed as

$$\begin{align} \begin{aligned} \omega_d = -i\mathbf{v_d}\cdot\nabla & = \underbrace{i\frac{c}{Z_\alpha}\left(\mu + \frac{m_\alpha v_{||}^2}{B_0}\right)\frac{1}{JB_0^2}\left[g\frac{\partial B_0}{\partial\theta}\frac{\partial}{\partial\psi}\right]}_\mathrm{\{I\}}\\ & \quad \underbrace{-n\frac{c}{Z_\alpha}\left(\mu + \frac{m_\alpha v_{||}^2}{B_0}\right)\frac{\partial B_0}{\partial\psi}}_\mathrm{\{II\}}\\ & \quad \underbrace{+ n\frac{cm_\alpha v_{||}^2}{Z_\alpha}\frac{1}{JB_0^2}\left(\frac{\partial g}{\partial\psi}q + \frac{\partial I}{\partial\psi}\right)}_\mathrm{\{III\}} \end{aligned}, \end{align} \tag{ 68 }$$

where term {I}, term {II} and term {III} represent the geodesic curvature, normal curvature and equilibrium current contributions respectively. In equation (68), we note that term {I} does not contribute to the precession frequency since the geodesic drift cancels over a bounce period, effect of term {III} is ignorable, while term {II} is equal to equation (67) in the limit of $\tau_b\to 0$, θ → 0 and $v_{||}\to 0$ (i.e. deeply-trapped limit)

$$\begin{align} \begin{aligned} \omega_{d0} = \overline{\omega}_{d,\tau_b\to 0} = - n\frac{c}{Z_\alpha}\left(\mu + \frac{m_\alpha v_{||}^2}{B_0}\right)\frac{\partial B_0}{\partial\psi}\Big\lvert_{\theta = 0} \end{aligned}, \end{align} \tag{ 69 }$$

which is used for the deeply-trapped approximation in equation (16), and only has ψ dependence, i.e. $\omega_{d0} = \omega_{d0}\left(\psi\right)$.

In order to elucidate the accuracy of $\omega_{d0}$, we compare equations (67) and (69) in both the analytic geometry [17] and EAST experimental geometry [10]. The safety factor q and magnetic shear $s = (r/q)(\mathrm{d}q/\mathrm{d}r)$ are given in figure 3. For analysis convenience, we use $\left(\psi_{bc},\theta_\mathrm{tip}\right)$ coordinates to represent trapped particles, where $\psi_{bc}\in[0,\psi_\mathrm{edge}]$ is the poloidal magnetic flux of banana orbit center and $\theta_\mathrm{tip} \in[0,2\pi]$ is the Boozer poloidal angle at the banana tip location (i.e. the mirror throat). Using the magnetic field strength at $\left(\psi_{bc},\theta_\mathrm{tip}\right)$ location, i.e. $B_{0,\mathrm{tip}} = B_0\left(\psi_{bc},\theta_\mathrm{tip}\right)$, we can express the pitch angle of trapped particle as $\lambda = B_a/B_{0,\mathrm{tip}}$, which varies in the range of $B_a/B_\mathrm{max}\lt\lambda\lt B_a/B_\mathrm{min}$. Substituting equation (65) into equation (67), the exact precession frequency $\overline{\omega}_d$ can be calculated for all trapped particles in terms of $\left(\psi_{bc},\theta_\mathrm{tip}\right)$. The ratio between $\overline{\omega}_d\left(\psi_{bc},\theta_\mathrm{tip}\right)$ and $\omega_{d0}\left(\psi_{bc}\right)$ in the analytic and EAST experimental geometries are shown in figures 4(a) and 5(a) respectively. First, it is seen that the numerical calculation of equation (67) agrees with the measurement from the realistic orbits by equations (61)–(64). Second, most areas on the low field sides are characterized with $\overline{\omega}_d/\omega_{d0}\sim 1$, thus $\omega_{d0}$ in equation (69) can well approximate the precession frequency for most trapped particles not only in the analytic geometry with low magnetic shear but also in the experimental geometry with broader range of magnetic shear, where the significance of magnetic shear term in equation (67) is decreased by the elongation effect.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, since the trapped particle fraction f_t determines the EE drive intensity through equations (21) and (22), one needs to adjust $\lambda_\mathrm{low}$ in equation (19) so that f_t only incorporates the contribution of trapped particles that satisfy $\overline{\omega}_d/\omega_{d0}\sim 1$, while the barely-trapped particles beyond the model capability should be excluded. Besides the constraint arising from the approximation on precession frequency, we note that the bounce average operations on electromagnetic fields in equation (8) are removed for integrating the EE moments, and corresponding validity regime of $|\theta|\ll1/|nq-m|$ becomes another constraint for $\lambda_\mathrm{low}$. Figures 4(b) and 5(b) show the 2D profile of f_t using $\lambda_\mathrm{low} = B_a/B_\mathrm{max}$ in the poloidal plane, which take into account all trapped particles and overestimate the deeply-trapped EE drive intensity. We then apply $\lambda_\mathrm{low} = 1$ in equation (19) to calculate f_t as shown in figures 4(c) and 5(c), which corresponds to the trapped particles with entire banana orbits on the low field side, in consistency with the validity regime of $\overline{\omega}_d/\omega_{d0}\sim 1$ in figures 4(a) and 5(a). Thus, the 2D function f_t calculated by using $\lambda_\mathrm{low} = 1$ in equation (19) correctly reflects the deeply-trapped EE drive intensity for modes peak at $k_{||}\sim 0$ and is applied in the following e-BAE simulations. For modes peak at moderate $k_{||}\unicode{x2A7D} 1/qR_0$ such as e-TAE, one needs to choose $\lambda_\mathrm{low}\lt1$ in equation (19) to calculate f_t in the more deeply-trapped regime (i.e. reduce θ_b ). For comparison, the widely used simple 1D form of $f_t \approx \sqrt{2\epsilon}$ assumes trapped particles distribute uniformly along the poloidal direction and might amplify the drive of deeply-trapped fraction. In addition, the omitted barely-trapped EE drive is in the higher order due to its small population.

To verify the global fluid-kinetic hybrid simulation model and corresponding numerical implementation in MAS code, we carry out simulations of e-BAE in analytic geometry based on a well-established benchmark case [17], where the safety factor profile and concentric-circular geometry are shown in figures 3(a) and 4 respectively. Protons are used for thermal ions with $Z_i = e$. The thermal electron, thermal ion and EE temperatures are uniform with $T_{i0} = T_{e0} = 500$ eV and $T_{h0} = 25 T_{e0}$, and the thermal electron density is uniform with $n_{e0} = 1.3\times 10^{14}$ cm⁻³. The EE density profile is described by $n_{h0} = 0.05n_{e0}\left[1+0.2\left(\mathrm{tanh}\left(0.26-\hat{\psi}\right)/0.06\right) -1.0\right]$, where $\hat{\psi} = \psi/\psi_w$ is the poloidal magnetic flux normalized by the wall value, and the reciprocal of EE density scale length $R_0/L_{nh} = 12.7$ peaks at q = 2 rational surface, where $L_{n,h} = \left(\nabla n_{h0}/n_{h0}\right)^{-1}$. The thermal ion density is then determined by the quasi-neutrality condition $Z_in_{i0} + q_e\left(n_{e0} + n_{h0}\right) = 0$. The on-axis magnetic field strength is $B_0 = 1.91$ T, the major radius is $R_0 = 0.65$ m, the minor radius is $a = 0.333R_0$, and the safety factor profile in figure 3(a) is $q = 1.797+0.8\hat{\psi}- 0.2\hat{\psi}^2$. Note that $\beta_{h}\sim \beta_{e}$ with these parameters, which is similar to fast ion pressure ordering for EI-driven AEs as shown in figure 3 of [34]. Nevertheless, here the EE density and temperature profiles are chosen for the dedicated benchmark study on e-BAE physics, and the EE-β threshold for AE excitation can be given by simulation and theory, while it is still difficult to identify the exact EE distribution in experiments.

The global mode structure and dispersion relation of n = 3 e-BAE from MAS simulations using $f_t\left(\lambda_\mathrm{low} = 1\right)$ are analyzed as follows. The 2D poloidal mode structure of electrostatic potential δφ is shown in figure 6(a1), which is characterized with a 'boomerang' shape. The dominant principal poloidal harmonic of m = 6 peaks at q = 2 rational surface, of which amplitude is much larger than the neighboring sideband harmonics of m = 5 and m = 7 that corresponds to a weakly ballooning structure as shown in figure 6(a2). Different from δφ, the poloidal mode structure of parallel vector potential $\delta A_{||}$ exhibits anti-ballooning feature as shown in figure 6(b1), where the real parts of m = 5 and m = 7 poloidal sidebands are comparable to the resonant m = 6 harmonic but with phase difference as shown in figure 6(b2). Moreover, the mode structure of $\Delta\phi = \delta\phi - \delta\psi$ is shown in figures 6(c1) and (c2), which represents the derivation from the ideal-MHD limit and reflects the polarization of each poloidal harmonics. The m = 6 harmonic is predominantly Alfvenic according to the small $\Delta\phi$ since the electrostatic component δφ and electromagnetic component δψ nearly cancel out, while $\Delta\phi$ perturbation is large in m = 5 and m = 7 sidebands which indicates the acoustic component becomes more important.

Zoom In Zoom Out Reset image size

It should be pointed out that there are ordering differences on constructing analytic equilibria between MAS and GTC early simulations in [17] though the parameter inputs are the same, and the EE mass $(m_{\mathrm{EE}})$ effect on e-BAE excitation is absent in [17] and awaits clarification. To investigate e-BAE physics in a more comprehensive way and careful benchmark MAS results, we carry out GTC simulations using current code version of GTC-4.6 with analytic equilibrium in the lowest order of inverse aspect ratio $\epsilon = r/R_0$ described in appendix E of [35] that are identical to MAS simulations in figure 6. GTC results of e-BAE mode structure and wave-particle resonance are shown in figure 7 with scaling of EE mass $m_{\mathrm{EE}}/m_i$, it is found that as $m_{\mathrm{EE}}/m_i$ decreases from $1/1$ to realistic $1/1836$, the δφ 2D mode structure exhibits a triangularity shape gradually, and the radial width of dominant m = 6 harmonic becomes narrower indicated by the grey shaded region in figures 7(a2)–(d2). In GTC simulations, the wave-particle resonance strength can be inferred from the intensity of EE entropy $\delta f_h^2$ in $\left(E,\lambda\right)$ phase space as shown in figures 7(a3)–(d3), it is seen that the dominant EE resonance moves from bounce-precession near λ = 1 to precession near $\lambda = 1+r/R_0\approx 1.16$ as $m_{\mathrm{EE}}$ value decreases to the realistic electron mass, which is consistent with recent theory that the bounce-averaged dynamics dominate EE resonances with low frequency MHD modes [13]. The corresponding e-BAE real frequency ω_r and growth rate γ at different $m_{\mathrm{EE}}/m_i$ values are shown in figure 8, it is noted that the variation of ω_r is small which is mostly determined by bulk plasmas, however, γ of $m_{\mathrm{EE}}/m_i = 1$ case significantly deviates from other cases using smaller $m_{\mathrm{EE}}/m_i$ values, because the resonance region in figure 7(a3) covers not only the precessional drift resonance but also the bounce-precession and transit resonances due to the unphysically large EE mass, while the resonance regions in figures 7(b3), (c3) and (d3) are dominated by precessional drift resonance and corresponding γ values are close with each other since trapped EE bounce-averaged dynamics only rely on temperature rather than mass [13]. Comparing the 2D mode structures and the 1D radial profiles of e-BAE m-harmonics between MAS simulation in figures 6(a1) and (a2) and GTC simulation in figures 7(d1) and (d2) using realistic $m_{\mathrm{EE}}/m_{i} = 1/1836$, it is seen that MAS eigenvalue result mostly agrees with GTC PIC simulation on δφ 2D mode structure and corresponding radial profile of each m-harmonic. However, MAS gives relatively lower amplitudes of m = 5 and m = 7 sideband harmonics than GTC which reflects EE drive strength, because MAS applies deeply-trapped approximation for calculating EE drive in figure 6 using trapped fraction on the side of $\lambda \gt \lambda_\mathrm{low} = 1$, while GTC simulation shows that a small part of $\delta f_h^2$ resonance structure in phase space extends to the barely-trapped region on the side of λ < 1 in figure 7(d3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regarding to the computational costs of MAS and GTC on e-BAE case above, MAS takes about 10 s with R = 100 radial grids and M = 3 coupled m-harmonics, and GTC runs 60 000 time steps by costing 55 GPU node hours on Perlmutter supercomputer (or 4000 CPU node hours) with mpsi = 128, mtheta = 512 and mtoroidal = 32 grids in the radial, poloidal and parallel directions and 75 particles per grid-cell by summing all species. The huge computational cost of GTC is caused by the smallness of realistic EE mass, resulting in the stringent restriction on time step by EE parallel transit time (i.e. $k_{||}v_{th,e}\Delta t\unicode{x2A7D} 1$) for initial value simulations [21]. The upgraded MAS code greatly improves the computational efficiency for modeling EE effects on low frequency MHD modes with $\omega\ll\Omega_{ci}$, which is useful for analyzing experimental observations and preparing data for training surrogate models in machine learning [36].

Next, we investigate parametric dependencies of e-BAE ω_r and γ on EE drive, namely, compare MAS and GTC results by varying $n_{h0}$ and $T_{h0}$ in amplitudes while remaining profiles unchanged. Considering the EE $\delta f_h^2$ resonance region in figure 7(d3) from GTC simulation, we perform MAS simulations with both $\lambda_{\mathrm{low}} = 1$ and $\lambda_{\mathrm{low}} = B_a/B_{\mathrm{max}}$ as the lower limits of pitch angle for trapped fraction f_t calculations using equation (19). From figures 9(a) and (b), it is seen that ω_r decreases and γ increases as $n_{h0}/n_{e0}$ varies from 0 to 0.1 in both MAS and GTC simulations. Specifically, MAS results using $\lambda_{\mathrm{low}} = 1$ agree with GTC results near the marginal stable regime of small $n_{h0}/n_{e0}$, and MAS results using $\lambda_{\mathrm{low}} = B_a/B_{\mathrm{max}}$ get close to GTC results as $n_{h0}/n_{e0}$ increases. The reason can be explained by comparing GTC $\delta f_h^2$ resonance regions in phase space between figures 9(c)–(e). For the weak EE drive case of $n_{h0}/n_{e0} = 0.03$, most intensity of $\delta f_h^2$ distributes on the side of λ > 1 where the deeply-trapped EE model using $\lambda_{\mathrm{low}} = 1$ can well approximate EE drive. However, the GTC $\delta f_h^2$ resonance regions gradually extend to λ < 1 side as $n_{h0}/n_{e0}$ increases, and the realistic EE drive is then underestimated in MAS using $\lambda_{\mathrm{low}} = 1$. On the other hand, MAS simulations using $\lambda_{\mathrm{low}} = B_a/B_{\mathrm{max}}$ overestimate EE drive that give higher γ than GTC because all trapped EEs are treated as deeply-trapped ones with precessional drift frequency that closely matches with the mode frequency. Similar conclusions are obtained for $T_{h0}/T_{e0}$ scan as shown in figure 10. It should be pointed out that the ellipticity-induced AE (EAE) can be excited by EEs and grows faster than e-BAE in GTC initial value simulations when $T_{h0}/T_{e0}\unicode{x2A7E} 40$, thus we cannot measure corresponding GTC results of e-BAE ω_r and γ for cases with $T_{h0}/T_{e0}\unicode{x2A7E} 40$ in figure 10(b). It is also seen that the GTC results of γ exhibit a transition from MAS results using $\lambda_{\mathrm{low}} = 1$ to $\lambda_{\mathrm{low}} = B_a/B_{\mathrm{max}}$ as $T_{h0}/T_{e0}$ increases in figure 10(b). The intersections of the black dashed lines and growth rate curves in figures 9 and 10(a), (b) indicate the excitation thresholds that EE drive just overcomes the continuum damping [37, 38] and radiative damping/Landau damping from bulk plasma [39–41], and it is clear that MAS model using $\lambda_{\mathrm{low}} = 1$ is sufficient for estimating these excitation thresholds, where the deeply-trapped EE play a dominant role in the marginal stable regime. In a short summary, MAS model considers the precessional drift resonance as the dominant EE resonance with AEs, which is appropriate and confirmed by GTC PIC simulations using e-BAE as an example. It is also worthy mentioning that these EE-driven AE processes have similarities with alpha particle driven AEs in future fusion reactor, of which EP orbits and wave–particle resonances are numerically investigated and compared in appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The symmetry breaking of AE mode structure due to non-ideal MHD effects, such as kinetic EPs induced anti-Hermitian part of dielectric tensor [42], have been widely observed in simulations [43, 44] and experiments [45, 46], which are characterized with twisted mode structures. Since the distortion of AE mode structure by EEs is much less explored compared to EIs in experiments, we also perform GTC first-principle simulations of i-BAEs in appendix B to check the symmetry property between i-BAE and e-BAE. By comparing the mode structures of i-BAE in figure 19 and e-BAE in figure 7, the main difference is the tip orientation of 'boomerang' (triangularity) shape, and it can be speculated that the underlying physics of EE non-perturbative effects on BAE mode structure is similar to EIs except for negative charge sign and small FOWs. To understand the role of EEs on 'boomerang' (triangularity) shape mode structure in figure 6(a) that is no longer up-down symmetric, four MAS simulation cases are performed with different EE physics as described in table 1. Case {I} corresponds to the Landau-fluid simulation of bulk plasma without any EE effects, and the poloidal mode structure is relatively up-down symmetric as shown in figure 11(a1), where the small distortion is due to the kinetic effects of bulk plasmas that lead to anti-Hermitian contribution. In case {II}, the EE-KPC term is added on top of case {I}, and the 2D mode structure is radially broadened and drastically twisted with obvious tails on both sides of mode rational surface in figure 11(b1), which indicates that the anti-Hermitian contribution from EE-KPC is much larger than bulk plasmas. For comparison, figure 11(c1) shows the 2D mode structure of case {III} that includes EE-IC term on top of case {I} instead of EE-KPC term, and it can be seen that EE-IC term has little impact on the symmetry breaking which is different from case {II}, but EE-IC term can broaden the radial width to a certain extent. In case {IV}, the non-perturbative effects of both EE-IC and EE-KPC terms are considered (figure 11(d1)), it is seen that the degree of symmetry breaking in case {IV} is between case {II} and case {III}, and the radial width is close to both case {II} and case {III}. Thus, there are two EE non-perturbative effects on e-BAE mode structure: (i) the EE-KPC term provides the dominant kinetic effects and induces the large ant-Hermitian part of dispersion relation, which is responsible for the distortion e-BAE mode structure; (ii) the EE-IC term represents the fluid-like convective response and enhances the Hermitian part, which not only compensates the symmetry breaking induced by EE-KPC term but also broadens the radial width through MHD interchange effect.

Zoom In Zoom Out Reset image size

Table 1. Four simulation cases of n = 3 e-BAE with different EE physics. Note that EE-IC and EE-KPC terms are given in equation (24), which represent the fluid and kinetic EE responses respectively.

Case	EE-IC	EE-KPC	$\omega_r\left(V_{Ap}/R_0\right)$	$\gamma\left(V_{Ap}/R_0\right)$
(I)	No	No	0.160	−0.007 07
(II)	No	Yes	0.175	0.004 96
(III)	Yes	No	0.134	−0.009 09
(IV)	Yes	Yes	0.149	0.009 04

The radial profile of phase angle θ_r can also reflect the mode structure symmetry, which has received interest from recent experiments [45, 46]. Here we combine θ_r radial profiles of cases {I}-{IV} with corresponding 2D mode structures to illustrate their underlying connections. Since the poloidal harmonics of e-BAE are weakly coupled, we choose the dominant m = 6 harmonic of δφ and calculate θ_r according to $\delta\phi_{m = 6} = |\delta\phi_{m = 6}|\mathrm{exp}(i\theta_r)$. In figure 12, the blue solid line represents θ_r , the red solid and red dashed lines represent the real and imaginary parts of $\delta\phi_{m = 6}$ respectively, and the radial domain of e-BAE eigenfunction with finite amplitude is marked as a gray shaded region. As shown in figures 12(a) and (c), θ_r profiles almost remain constant in the e-BAE region of cases {I} and {III}, which correspond to the relatively symmetric mode structures in figures 11(a1) and (c1) respectively. In contrast, θ_r profiles show large variations in the e-BAE region of cases {II} and {IV} due to the anti-Hermitian EE-KPC term that breaks the poloidal up-down symmetry. However, the poloidal phase shifts at different radial locations are approximately symmetric with respect to the e-BAE peak as shown in figures 12(b) and (d) (because the EE gradient profile peaks at the mode rational surface in our simulations), which still indicate the radial symmetry of 'boomerang' shape mode structures characterized with two symmetric tails as shown in figures 11(b1) and (d1). In addition, we show the real frequencies, radial positions and mode widths of cases {I}–{IV} on corresponding n = 3 continuous spectra in figure 13, where the Alfvenic and acoustic continua are calculated in MAS code based on an ideal full-MHD model as described in appendix B of [25].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Further, in order to clarify the competition between EE excitation and bulk plasma damping, we compute the EE drive $\gamma_{\mathrm{drive}}$ in the absence of dissipation by dropping Landau damping terms associated to $|k_{||}|$ in equations (25)–(27), and compute the bulk plasma damping rate $\gamma_{\mathrm{damp}}$ by isolating EE drive effect (i.e. only keeping the real part of EE response functions in equations (21) and (22) while dropping the imaginary part), and $\gamma_{\mathrm{drive}}$ and $\gamma_{\mathrm{damp}}$ dependencies on $n_{h0}$ and $T_{h0}$ are shown in figure 14 for EE trapped fraction f_t calculated with both $\lambda_{\mathrm{low}} = 1$ and $\lambda_{\mathrm{low}} = B_a/B_{\mathrm{max}}$. It should be noted that both bulk plasma and EE non-perturbative effects, namely, modifications on mode structure and real frequency, are kept for calculating $\gamma_{\mathrm{drive}}$ and/or $\gamma_{\mathrm{damp}}$. For example, $\gamma_{\mathrm{damp}}$ shows relatively strong dependency on $T_{h0}$ in figure 14(b), which indicates that EE non-perturbative effect can quantitatively affect the bulk plasma damping on e-BAE through altering the mode structure and real frequency, and thus becomes necessary for accurate damping and growth rate calculation. Readers may also notice that the change of $\gamma_{\mathrm{damp}}$ by varying $T_{h0}$ becomes weaker in figure 14(d), the reason can be understood as follows. We recall that the EE-IC term incorporates both passing and trapped EE contributions (i.e., adiabatic fluid convection) in MAS model, while $f_t$ only acts on EE-KPC term in equation (24), and the non-perturbative corrections on real frequency are opposite in sign between EE-IC term and EE-KPC term as illustrated in figure 13, thus larger $f_t$ using $\lambda_{low} = B_a/B_{max}$ in figure 14(d) gives rise to a larger cancellation of EE non-perturbative effect on real frequency, and the resonance condition for damping is less altered. Meanwhile, the $\gamma_{\mathrm{drive}} - |\gamma_{\mathrm{damp}}|$ agrees well with the gross growth rate γ using comprehensive model with both EE drive and bulk plasma damping, which confirms the correctness of our method for obtaining $\gamma_{\mathrm{drive}}$ and $\gamma_{\mathrm{damp}}$ that covers not only the unstable regime of $|\gamma_{\mathrm{damp}}|\ll\gamma_{\mathrm{drive}}$ but also the marginally stable regime of $|\gamma_{\mathrm{damp}}|\sim\gamma_{\mathrm{drive}}$, and thus provides useful damping and drive information for linear stability and further nonlinear dynamics analyses [47].

Zoom In Zoom Out Reset image size

Finally, we present MAS simulation of e-BAEs in experimental geometry based on EAST discharge #82589 at 4090 ms, of which q profile and magnetic geometry are shown in figures 3(b) and 5. Within the ECRH and LHCD in EAST tokamak, multiple AEs with toroidal mode numbers n = 1 and n = 2 propagating along electron diamagnetic drift direction are observed in a wide frequency range of $[20\,\textrm{kHz}, 300\,\textrm{kHz}]$. Meanwhile, the EE population is enhanced during these AE activities based on the diagnosis of hard x-ray photon counts. Since EI effect is ignorable without applying ion cyclotron resonance heating (ICRH) and neutral-beam injection (NBI), the AEs in this EAST discharge are destabilized by EEs as reported in [10]. The thermal plasma profiles and corresponding Alfven continua are shown in figures 6 and 7 of [10], here we focus on $\left(m/n\right) = \left(4/1\right)$ BAE mode and analyze its stability in $n_{h0}-T_{h0}$ space, the EE density profile is nonuniform given by $n_{h0} = C\times n_{e0,a}\left[1+0.3\left(\mathrm{tanh}\left(0.85-\hat{\psi}\right)/0.06\right) -1.0\right]$, where $C\ll 1$ is a constant and $n_{e0,a}$ represents the thermal electron density at magnetic axis, and the EE temperature profile $T_{h0}$ is uniform. The MAS simulations are performed with $\lambda_{\mathrm{low}} = 1$, and the parametric dependencies of e-BAE ω_r and γ on $n_{h0}$ and $T_{h0}$ are shown in figure 15. It is seen that as $\beta_h = 8\pi n_{h0}T_{h0}/B_0^2$ increases, ω_r slightly decreases from the continuum accumulation point (CAP) near 32 kHz, while $\gamma/\omega_r$ quickly increases from negative damping to positive growth, where the cyan solid line indicates the marginal stable regime in figure 15(b). The $T_{h0}/T_{e0,a}$ range for e-BAE excitation is $[15, 35]$ along the cyan solid line that corresponds to [63 keV, 149 keV], which partially overlaps with experimental measurements in figure 3 of [10], and there exists a minimal value of $\beta_{h,\mathrm{crit}} \approx 0.76$ $\beta_{e0,a}$ (where $\beta_{e0,a}$ is on-axis thermal electron pressure ratio) for e-BAE excitation indicated by the red pentagram in figure 15(b). Given the fact that precession frequency is proportional to both n number and energy, both upper and lower limits of EE temperature range for $n=2$ e-BAE excitation decrease by a half than that of $n=1$ case, which agrees with the experimental measurements of $T_{h0} < 63$ keV part. The $m/n = 4/1$ BAE mode structure in the absence of EE drive is shown in figure 16, and the $m/n = 4/1$ e-BAE mode structure with EE drive at minimal $\beta_{h,\mathrm{crit}}$ location in figure 15(b) is shown in figure 17. The e-BAE exhibits larger amplitudes of m = 3 and m = 5 sideband harmonics and radially broader mode structure, compared to BAE in the absence of EE drive. Note that e-BAE in figure 17 is more radially localized compared to e-BAE in figure 6, because the magnetic shear in EAST is large and the neighboring rational surfaces are more close to each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size