# Comprehensive linear model for the n = m = 1 fishbone kinetic-MHD instability

G. Brochard, R. Dümont, X. Garbet, H. Lutjens, T. Nicolas, F. Orain

► To cite this version:

G. Brochard, R. Dümont, X. Garbet, H. Lutjens, T. Nicolas, et al. Comprehensive linear model for the n = m = 1 fishbone kinetic-MHD instability. Journal of Physics: Conference Series, 2018, 1125, pp.012003. 10.1088/1742-6596/1125/1/012003. [hal-03319164](https://hal.science/hal-03319164)

HAL Id: hal-03319164

https://hal.science/hal-03319164

Submitted on 10 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comprehensive linear model for the $n = m = 1$ fishbone kinetic-MHD instability

G. Brochard<fnref n="1" />, R. Dumont<fnref n="1" />, X. Garbet<fnref n="1" />, H. Lütjens<fnref n="2" />, T. Nicolas<fnref n="2" />, F. Orain<fnref n="2" />

E-mail: guillaume.brochard@cea.fr

**Abstract.** A comprehensive linear model for the interaction between the MHD $n = m = 1$ internal kink and fast particles, based on the Hamiltonian angle-action formalism, is derived. On the basis of this model, a linear code, MHD-K, that solves partially analytically and numerically the kinetic-MHD dispersion relation non-perturbatively is presented. The impact of passing fast particles on the fishbone is shown to be an essential drive of the instability, where previous models highlighted only trapped particles as the driver of the kinetic-MHD instability. Resonant planes in phase space are presented, showing multiple resonant branches for both trapped and passing particles.

## 1. Introduction

The kinetic-MHD fishbone instability has been studied extensively for the past three decades, from the experimental, theoretical and numerical points of views. First results of MHD bursts with the typical fishbone-like signature on magnetic oscillations were found at the PDX tokamak [1]. Since then, they have been observed in many other experiments [2][3] for different phase space distributions of energetic particles. The relevance of studying this instability rests upon the fast particles transport, induced by the fishbone instability. The fast particles, such as the fusion born alpha particles, are transported out of the plasma core. This decreases the plasma performances due to a partial loss of the alpha heating, mandatory for steady-state plasma operations. Furthermore, such a transport can trigger disruptives MHD instabilities, such as NTMs (Neoclassical Tearing Modes). Due to the fact that future fusion devices such as ITER and DEMO are expected to produce a large quantity of alpha particles, understanding the interplay between kinetic suprathermic population and MHD both in its linear and non-linear phase is crucial for conducting successful experiments. While the linear phase can be investigated analytically for simple models, or with little numerical cost for more complete models, the study of the non-linear phase of the fishbone instability requires the use of complex numerical codes such as MEGA [4], M3D-K [5][6], XTOR-K[7][8]. However, for realistic non-linear predictions, such codes need to be firstly validated against complete linear model during their linear phase. The goal of this study is to provide analytical/numerical linear tools to perform such a validation, by deriving a complete linear model for the fishbone instability, and a linear code to solve the implicit dispersion relation. Using the formalism from [9], a more complete model is needed because, so far, theoreticallinear models [10],[11],[12],[13], primarily identified the precessionnal frequency of trapped energetic particles $\omega_d$, as the drive of a resonant mechanism with the $n = m = 1$ internal kink mode, inducing the fishbone instability. $v_d$ is the magnetic drift velocity. It is argued in this paper that passing fast particles need to be considered as an equivalent drive to the fishbone instability. In later models [14] the contribution of passing particles has been taken into account, however the impact on the dynamics of the fishbone has often been discarded, in some linear codes [5][15]. Moreover, some linear codes [16] compute perturbatively the fast particles contribution to the internal kink dispersion relation. Such an approximation holds for kinetic instabilities such as TAEs (Toroidal Alfven Eigenmodes). However, in the case of the fishbone, the growth rate of the instability can be up to one order of magnitude higher for a given magnetic equilibrium on the fishbone branch than it is on the kink branch without

<figure><img src="image_2.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France</footnote>

<footnote n="2">² CPhT, Ecole Polytechnique, F-91128 Palaiseau cedex, France</footnote>
fast particles. That is why a non-perturbative treatment of the kinetic dispersion relation, implicit in complex frequency as it will be shown, is needed. This paper is organized as follows, in a first part, a comprehensive linear model for the fishbone instability is derived, through a Hamiltonian angle-action formalism that will be presented as well. In a second part, the linear code MHD-K that solves partly analytically and numerically the fishbone dispersion relation is presented, as well as linear results showing the relevance of passing particles on the instability due to multiple resonant planes in phase space associated with both passing and trapped particles.

## 2. Derivation of the linear model

### 2.1. Kinetic dispersion relation of the internal kink

Most linear models have used the energetic principle [17] to work out an expression for the fishbone dispersion relation with the additional functional introduced by kinetic populations : $\delta W_K = d^3\mathbf{x}d^3\mathbf{v} m_h \xi \cdot \nabla \cdot (\mathbf{v} \otimes \mathbf{v})\tilde{f}_h$, $\tilde{f}_h$ the linear distribution response of energetic particles and $\xi$ the MHD displacement. This kinetic term comes from the extra stress tensor $\Pi_h$ due to fast particles in the momentum equation. Following what has been done in [18] and also in [19][20], the dispersion relation for the bi-fluid resistive internal kink mode $n = 1, m = 1$ is obtained by matching the expressions for $\xi_r$ obtained from the linearized MHD equations, inside and outside of the inertial layer $q = 1$, yielding

$$D(n_{h0}, \omega) = \delta I(\omega) - i\omega_A[\lambda_H + \lambda_K(n_{h0}, \omega)] = 0 \quad (1)$$

$\omega \in \mathbb{C}$, with the inertial contribution $\delta I$ given as

$$\frac{\delta I(\omega) = [3\omega(\omega - \omega_{*i})] I_R(\omega), \quad I_R(\omega) = \frac{8\Gamma \frac{\Lambda^{3/2} + 5}{4}}{\Lambda^{9/4} \Gamma \frac{\Lambda^{3/2} - 1}{4}}, \quad \Lambda = [\omega(\omega - \omega_{*i})(\omega - \omega_{*e})] \tau_A s_0^2}{\frac{1}{3} \quad \frac{S}{\Gamma} \quad \frac{1}{3}} \quad (2)$$

where $\omega_{*i}, \omega_{*e}$ stands for the ion/electron-bulk diamagnetic frequency, $\omega$ the complex MHD mode frequency, $S$ the Lundquist number, $s_0$ the magnetic shear at $q = 1$, $\omega_A$ the Alfvén pulsation, $n_{h0}$ the density of fast particles on the magnetic axis, $\lambda_H = \gamma_{MHD}\tau_A$ the contribution of the bulk plasma, $\tau_A$ being the Alfvén time, and $\gamma_{MHD}$ is the growth rate of the instability without fast particles, constant for a given MHD equilibrium. $\lambda_K$ is a normalization of $\delta W_K$ given by $\lambda_K = \mu_0\delta W_K / (2\pi\xi_0^2 s_0 R_0 B_0^2)$, with $R_0$ the major radius, $B_{\theta,0}$ the poloidal field at $q = 1$ and $\xi_0$ the constant radial MHD displacement. Here the eigenfunction $\xi$ has been assumed to be a simple step function, null outside $q = 1$. The way (1) is solved will be discussed in the second part of the article. Firstly, the derivation of the kinetic term $\delta W_K$ is presented. For this purpose, one needs to solve the Vlasov equation to work out an expression for the perturbed distribution of fast particles $\tilde{f}_h$.

### 2.2. Resolution of Vlasov equation

Let $F_h$ be the distribution function of fast particles such as $F_h = F_{eq,h} + \tilde{f}_h$, with $F_{eq,h}$ being the equilibrium distribution function. Vlasov's equation can be written as, assuming no collisions

$$d_t F = \partial_t F - \{H, F\} = 0 \quad (3)$$

where $\{\cdot, \cdot\}$ are the the Poisson's brackets defined in an angle/action formalism [21] as

$$\{A, B\} = \frac{\partial A}{\partial \alpha} \frac{\partial B}{\partial \dot{J}} \frac{\partial B}{\partial \dot{\alpha}} \frac{\partial A}{\partial J} \quad (4)$$

In a tokamak configuration, one can identify three invariants of motion $J_i$ related to a set of conjugate fast varying angles $\alpha_i$ [22],[23]. The conjugate coordinates are given as follows

* $(\alpha_1, J_1)$ is related to the cyclotron motion, where $\alpha_1$ stands for the gyroangle and $J_1$ is proportional to the magnetic moment $\mu$

* $(\alpha_2, J_2)$ is related to the poloidal motion, where $\alpha_2$ stands for the bounce angle and $J_2$ is a function of the two other invariants and the energy

* $(\alpha_3, J_3)$ is related to the toroidal motion, where $J_3 = P_\phi = mRv_\phi - e\psi(r)$ is the toroidal kinetic momentum, $m$ and $e$ the particle's mass and charge, $v_\phi$ its toroidal velocity and $\psi$ the poloidal magnetic flux.
This set of coordinates has been built so that the equilibrium quantities are only functions of **J**, therefore the linearized equations of motion regarding the equilibrium Hamiltonian $H_{eq}$, where $H(\alpha, \mathbf{J}) = H_{eq}(\mathbf{J}) + \tilde{h}(\alpha, \mathbf{J})$, $\tilde{h}$ the perturbed hamiltonian, simply read

$$ \frac{d\alpha}{dt} = \frac{\partial H_{eq}}{\partial \mathbf{J}} = \Omega, \qquad \frac{d\mathbf{J}}{dt} = -\frac{\partial H_{eq}}{\partial \alpha} = 0 \qquad (5) $$

Therefore, the unperturbed motion is characterized by three frequencies. $\Omega_1 = \omega_c$ is the gyrofrequency and $\Omega_2, \Omega_3$ are respectively the poloidal and toroidal transit frequency given by [21]

$$ \Omega_2 = \omega_b = 2\pi \int_{-\vartheta_0}^{\vartheta_0} \frac{d\vartheta}{i\mathbf{v}_{\|}} e^{-\frac{1}{2}\int_{-\vartheta_0}^{\vartheta_0} \frac{d\vartheta}{i\mathbf{v}_{\|}}} , \quad \Omega_3 = \delta_p q \omega_b + \omega_d \qquad (6) $$

where $\vartheta_0$ is the bounce angle, $l$ the Jacobian for toroidal coordinates $(\psi, \vartheta, \phi)$, $\vartheta$ the geometric poloidal angle, $v_{\|}$ the parallel velocity, $\delta_p$ is 1 for passing particles and 0 for trapped particles and $\omega_d$ the precessional frequency. Linearizing (3), one gets

$$ \partial_t \tilde{f}_h - \{\tilde{h}, F_{eq,h}\} - \{H_{eq}, \tilde{f}_h\} = 0 \qquad (7) $$

reducing, according to the definitions above, to

$$ \partial_t \tilde{f}_h - \frac{\partial \tilde{h}}{\partial \alpha} \cdot \frac{\partial F_{eq}}{\partial \mathbf{J}} + \frac{\partial \tilde{f}_h}{\partial \alpha} \cdot \frac{\partial H_{eq}}{\partial \mathbf{J}} = 0 \qquad (8) $$

since $F_{eq}$ does not depend on angle variables, as it is an equilibrium quantity. Performing a Fourier expansion on the perturbed quantities, according to the defined angle/action variables, one has

$$ \tilde{g} = \sum_{\mathbf{n}} g_{\mathbf{n}\omega} e^{i(\mathbf{n} \cdot \alpha - \omega t)} \qquad (9) $$

where $\mathbf{n} = (n_1, n_2, n_3)$ is linked to the angle/action formalism. Therefore, one readily obtains

$$ \tilde{f}_{\mathbf{n}\omega} = - \frac{\mathbf{n} \cdot \partial F_{eq,h} / \partial \mathbf{J}}{\omega - \mathbf{n} \cdot \Omega} \tilde{h}_{\mathbf{n}\omega} \qquad (10) $$

Considering that $F_{eq,h}$ is solely a function of $(H_{eq}, \mu, P_{\phi})$, its derivative with respect to **J** expands as

$$ \mathbf{n} \cdot \frac{\partial F_{eq,h}}{\partial \mathbf{J}} = \mathbf{n} \cdot \frac{\partial H_{eq}}{\partial \mathbf{J}} \frac{\partial F_{eq}}{\partial H_{eq}} + \frac{\partial \mu}{\partial \mathbf{J}} \frac{\partial F_{eq}}{\partial \mu} + \frac{\partial P_{\phi}}{\partial \mathbf{J}} \frac{\partial F_{eq}}{\partial P_{\phi}} = \frac{\partial F_{eq,h}}{\partial H_{eq}} \frac{\mathbf{h}}{\mathbf{n} \cdot \Omega - \omega} + \omega - \omega + \omega_{*} \quad i \qquad (11) $$

where $\omega$ is artificially introduced to later dissociate $\tilde{f}$ in resonant and fluid parts. $\omega_+$ and $\omega_*$ are defined as

$$ \omega = n \frac{e_s \partial F_{eq,h} / \partial \mu}{\omega} , \qquad \omega = -n \frac{\partial F_{eq,h} / \partial P_{\phi}}{\partial F_{eq,h} / \partial H_{eq}} \qquad (12) \\ + \quad m \partial F_{eq,h} / \partial H_{eq} \qquad * \quad 3 \partial F_{eq,h} / \partial H_{eq} $$

with $e_s$ the particle's charge and $m$ the particle's mass. Considering the perturbed Hamiltonian, derived from the electromagnetic Lagrangian, assuming that the perpendicular potential vector $A_{\perp}$ can be discarded

$$ \tilde{h} = e_s (\varphi - v_{\|} A_{\|}) \qquad (13) $$

with $\varphi$ the electric potential. Following the ideal MHD constraint on the system, $E_{\|} = 0$, and introducing the time integral of the electrical potential $\chi(\mathbf{x}, t) = \int_0^t \varphi(\mathbf{x}, t') dt'$, it leads to

$$ E_{\|} = - \nabla_{\|} \varphi - \partial_t A_{\|} = 0 \Leftrightarrow A_{\|} = - \nabla_{\|} \chi \qquad (14) $$

The perturbed Hamiltonian can then be expressed as

$$ \tilde{h} = e_s [\partial_t \chi + v_{\|} \nabla_{\|} \chi] = e_s [d_t \chi - (v_{\|} \cdot \nabla) \chi] \qquad (15) $$
where $d_i\chi = (\partial_t + \mathbf{v} \cdot \nabla)\chi, \mathbf{v} = \mathbf{v}_{\parallel}\mathbf{b} + \mathbf{v}_d, \mathbf{v}_d$ being the drift velocity of fast particles. $d_i\chi$ can be expressed, using the angle-action Hamiltonian approach, as

$$d_i\chi = \partial_t\chi - \{H_{eq}, \chi\} = \partial_t\chi + \Omega \cdot \frac{\partial\chi}{\partial\alpha} \qquad (16)$$

Therefore, plugging (16) into (15) and Fourier developing $\tilde{h}$, one gets

$$\tilde{h}_{n\omega} = i(n \cdot \Omega - \omega)e_s\chi_{n\omega} - [e_s(v_d \cdot \nabla)\chi]_{n\omega} \qquad (17)$$

Then, inserting (17),(11) into (10), one obtains $\tilde{f}_{n\omega} = \tilde{f}_{n\omega}^{res} + \tilde{f}_{n\omega}^{fl}$, with

$$\tilde{f}^{res} = \frac{\partial F_{eq}}{\partial H_{eq}} \frac{\omega - \omega_{+} - \omega_{*}}{[e_s(\mathbf{v}_d \cdot \nabla)\chi]} , \quad \tilde{f}^{fl} = \frac{\partial F_{eq}}{\partial H_{eq}} i(n \cdot \Omega - \omega - \omega)e\chi - [e(\mathbf{v} \cdot \nabla)\chi] \quad (18)$$

Let us now introduce the ordering between the different resonant, diamagnetic and MHD frequencies for energetic particles. $\Omega_1$ being generally several orders of magnitudes above $\Omega_2, \Omega_3$ and the pulsation of the internal kink $\omega$, fast particles gyrofrequency cannot react resonantly with the MHD mode. It is then possible to reduce the number of resonances considering $n_1 = 0$. It implies also that $\omega_+ = 0$. Given that perturbed quantities such as the electrostatic potential are assumed to have, due to the internal kink, the following form $\varphi(r, \vartheta, \varphi) = \varphi(r)e^{-i(\vartheta-\varphi)}$, their Fourier transform using angle-action formalism gives

$$\varphi_n = \frac{1}{(2\pi)^2} \int d\alpha_2 d\alpha_3 \varphi(r)e^{-i(\vartheta-\varphi+n_3\alpha_3+n_2\alpha_2)} \qquad (19)$$

Since $\alpha_3 = \phi - q\hat{\vartheta}$, with $\hat{\vartheta}$ a periodic function of $\alpha_2$ defined in [21], and that $\vartheta = \hat{\vartheta} + \delta_p\alpha_2$, the only non vanishing component of $\varphi_n$ following the integration over $\alpha_3$ is $n_3 = n = 1$. Therefore, the general resonance condition reads

$$\omega - \sum_{n_2} \sigma_{\parallel}(n_2 + \delta_p q n)\omega_b - n\omega_d = 0 \qquad (20)$$

$\sigma_{\parallel} = \pm 1$ for co/counter passing particles. Since the frequencies $\omega_b, \omega_d$ depend on the three invariants of motion $(E, P_{\phi}, \mu)$, the solutions of this resonance condition are surfaces in phase space. After numerical studies, which will be presented in part 3.3, for typical fast particles distribution and ITER-like MHD equilibrium, only three surfaces, one for trapped particles with $n_2 = 0$ and two for passing ones with $n_2 = -1$, mainly contribute to the resonant term, yielding for the $n = m = 1$ internal kink

$$\omega + \sigma_{\parallel}\delta_p(1 - q)\omega_b(E, P_{\phi}, \mu) - \omega_d(E, P_{\phi}, \mu) = 0 \qquad (21)$$

Regarding the fluid term, it is not necessary to use the angle/action formalism, this formalism being better adapted to the study of resonant terms. It is therefore more convenient to express it in real space. Noticing that in Fourier space $\mathbf{v} \cdot \nabla = i\mathbf{n} \cdot \Omega$ and that $\{\chi, P_{\phi}\} = i\eta_3\chi_{n,\omega}$

$$\tilde{f}^{fl} = e_s \frac{\partial F_{eq}}{\partial H_{eq}} \mathbf{v}_{\parallel} \nabla_{\parallel} \chi - \omega_{*} \{\chi, P_{\phi}\} / n_3 \qquad (22)$$

The first term of this expression is directly linked to the particle's tearing term in [9]. As a current term, it needs to be incorporated in $\delta W_{MHD}$, since it uses the total current, as shown in eq. (85-86) of [14].

In order to derive explicitly $\tilde{f}$, an expression of $\chi$ as a function of $\xi_0$, the order 0 MHD displacement, is needed. For $n = m = 1$, it reads $\xi = \xi_0 H(r_{q=1} - \tilde{r})e^{-i(\vartheta-\varphi+\omega t)}e_r$ at lowest order, where H is the Heaviside function and $\tilde{r}$ a radial position linked to a specific choice of referent flux surface $\tilde{\psi}$ [17][29]. Following reduced MHD assumption, at lowest order in $\epsilon = \tilde{r}/R$ and assuming a toroidal equilibrium with concentric circular surfaces without Shafranov shift

$$\mathbf{v}_{\perp} = \mathbf{v}_{E \times B} \Leftrightarrow -i\omega\xi_0 B_0 e^{-i(\vartheta-\varphi+\omega t)} = i\omega\partial_{\vartheta}\chi/\tilde{r} \Leftrightarrow \chi = -i\xi_0 B_0 \tilde{r} e^{-i(\vartheta-\varphi+\omega t)} \qquad (23)$$
It is possible to obtain explicit expressions for both the resonant and fluid perturbed distribution using the definition of the drift velocity

$$ \mathbf{v}_d = \frac{2E}{m\omega_c B} \mathbf{B} \times \left(1 - \frac{\lambda}{H} \mathbf{K} + \frac{\lambda}{2H} \nabla \ln B\right) \quad (24) $$

where $\lambda = \mu_0 B_0 / H_{eq}$ stands for the generalized pitch angle, $\mu_0$ being the order 0 magnetic moment and $B_0$ the magnetic field at $\bar{r} = 0$, $H = R / R_0 = 1 + \bar{r} / R_0 \cos \vartheta$, $\omega_c$ the particle's gyrofrequency and $\mathbf{K}$ the magnetic curvature. At low $\vartheta$, this expression reduces to

$$ \mathbf{v}_d = \frac{E\sigma}{m\Omega B H} \mathbf{B} \times \nabla \ln B \quad (25) $$

with $B = B_0 R_0 / R$, $\sigma = 2 - \lambda / H$ and $\mathbf{B} = (B_0 / H)(\mathbf{e}_\varphi + (\bar{r} / R q)\mathbf{e}_\vartheta)$. Therefore

$$ \nabla \ln B = \frac{h \sin \vartheta}{R_0} \mathbf{e}_\vartheta - \frac{\cos \vartheta}{R_0} \mathbf{e}_r \quad (26) $$

Then

$$ \mathbf{v}_d \cdot \nabla = - \frac{\sigma E}{e_s B_0 R_0} \begin{array}{c} \sin \vartheta \partial_r \\ \cos \vartheta \partial_\vartheta \end{array} \quad (27) $$

For the resonant term, one has

$$ (\mathbf{v}_d \cdot \nabla)\chi = \frac{\sigma E \xi_0}{e_s R_0} e^{i(\alpha_3 + q\hat{\vartheta} - \omega t)} \quad (28) $$

The Fourier coefficient being given by the inverse Fourier transformation

$$ g_{n\omega} = \int \frac{d\alpha_1}{2\pi} \frac{d\alpha_2}{2\pi} \frac{d\alpha_3}{2\pi} e^{-i\mathbf{n} \cdot \mathbf{\alpha}} g \quad (29) $$

two cases need to be considered in computing the resonant term. Firstly, for trapped particles, as will be shown numerically in part 3.3 for the n=m=1 internal kink, only the triplet $\mathbf{n}_t = (0, 0, 1)$ contributes, therefore

$$ [e_s (\mathbf{v}_d \cdot \nabla) \chi]_{\mathbf{n}_t \omega} = \frac{E \xi_0}{R} \int \frac{d\alpha_2}{2\pi} \sigma e^{iq\vartheta} = \frac{E \xi_0}{R} \sigma I^q \quad (30) $$

where $I_q = \langle \cos \vartheta \rangle_{\alpha_2} + (1 - q) \langle \vartheta \sin \vartheta \rangle_{\alpha_2}$, whose derivation is given in annex. For passing particles, only the triplet $\mathbf{n}_p = (0, -1, 1)$ contributes, therefore

$$ [e_s (\mathbf{v}_d \cdot \nabla) \chi]_{\mathbf{n}_p \omega} = \frac{E \xi_0}{R} \int \frac{d\alpha_2}{2\pi} \sigma e^{iq\vartheta + (1-q)\alpha_2} = \frac{E \xi_0}{R} \sigma I_{q,p} \quad (31) $$

Given the $\alpha_2$ dependency in $I_{q,p}$, the obtention of an analytical expression for it is quite intricated. However, performing a serie expansion of $\exp[(1-q)\alpha_2]$, where $1-q$ is considered as a small parameter given that the studied MHD instability is the internal kink, trapped and passing contributions are identical at lowest order. Therefore the analytical expression for $I_q$ is used for both trapped and passing particles. The resonant perturbed distribution function then reads

$$ \tilde{f}_{n\omega}^{res} = E \frac{\partial F_{eq}}{\partial \omega} \quad \omega - \omega^* \quad \frac{\xi_0}{I} (r, \lambda) \quad (32) \\ \partial H_{eq} \omega + \delta_p \sigma_{\parallel} (1-q) \omega_b - \omega_d R_0^q $$

Regarding the fluid term, using the expression of Poisson's brackets in regular toroidal coordinates [24], one finds

$$ \{\chi, P_{\phi}\} = - \frac{\mathbf{b}}{eB_0} \cdot \nabla \chi \times \nabla P_{\phi} = - \frac{\bar{r}}{q} \cdot (\mathbf{b} \times \nabla \chi) \quad (33) $$

where $P_{\phi} = \text{cst} = -e\bar{\psi} = -e\bar{r}^2 B_0 / 2q$, assuming thin orbits around the magnetic surfaces and $q(\bar{r})$ almost constant inside $\bar{r} = 1$. The choice $\bar{\psi} = \psi + mRv_{\phi}/Ze$ is explained in the next section. The fluid contribution, or interchange contribution by comparison with [9], yields
$$ \tilde{f}^{\text{int}} = -\xi_0 \frac{\partial F_{\text{eq}}}{\partial \bar{r}} e^{i(\phi - \vartheta)} \qquad (34) $$

IOP Conf. Series: Journal of Physics: Conf. Series **1125** (2018), 012003, doi:10.1088/1742-6596/1125/1/012003 where one can see that the adiabatic contribution from [12] is retrieved, as well as the more general expression in eq. (71) of [14] in the limit of thin particles orbits.
### 2.3. Derivation of the toroidal drift frequency $\omega_d$

To retrieve an expression for $\omega_d$, one needs to express the toroidal angle as $\phi = \omega_3 t + F(\alpha_2)$, where $\omega_3$ is the pulsation associated with the third angle of the angle-action formalism, and $F$ is a periodic function of $\alpha_2$. Given the angle $\alpha$ from the Clebsch magnetic field expression, $\alpha = \phi - q(\psi)\vartheta$; and that any angle $\xi$ respects $\xi = \mathbf{v} \cdot \nabla = (v_\| \nabla_\| + v_d \cdot \nabla)\xi$

$$ \dot{\alpha} = \mathbf{v}_d \cdot \nabla \phi - q(\bar{\psi})\vartheta - \frac{dq}{d\psi}(\bar{\psi})\hat{\psi}\vartheta \quad (35) $$

where $\bar{\psi}$ stands for the referent flux surface for a given particle orbit, an invariant as a function of invariants. $\hat{\psi}$ stands for the excursion from $\bar{\psi}$ along the orbit, such as $\psi = \bar{\psi} + \hat{\psi}$. Given that $\bar{\psi}$ is an unperturbed quantity, and that $\hat{\psi}$ is a first order quantity, as well as $\mathbf{v}_d \cdot \nabla\vartheta$, at leading order

$$ \dot{\alpha} = \mathbf{v}_d \cdot \nabla \phi - q(\bar{\psi})\mathbf{v}_d \cdot \nabla \vartheta - \frac{dq}{d\psi}(\bar{\psi})\vartheta \mathbf{v} \cdot \nabla \hat{\psi} \quad (36) $$

Noticing that $\vartheta \mathbf{v}_d \cdot \nabla \hat{\psi} = \vartheta(d\hat{\psi}/dt) = d(\vartheta \hat{\psi})/dt - (d\vartheta/dt)\hat{\psi}$

$$ \dot{\alpha} = \mathbf{v}_d \cdot \nabla \phi - q(\bar{\psi})\mathbf{v} \cdot \nabla \vartheta + \frac{dq}{d\psi}(\bar{\psi})\hat{\psi} \frac{d\vartheta}{dt} - \frac{dq}{d\psi}(\bar{\psi}) \frac{d(\hat{\psi}\vartheta)}{dt} \quad (37) $$

Then, performing an averaging over $\alpha_2$, and a series expansion in $q(\psi)$ regarding the toroidal angle, $\phi$ can be written as

$$ \phi = (\omega_d + \delta_p q(\bar{\psi})\omega_b)t + \delta_p q(\bar{\psi})\hat{\psi}(\alpha_2) + \chi(\alpha_2) \quad (38) $$

where the toroidal drift frequency is given by

$$ \omega_d = \mathbf{v}_d \cdot \nabla \phi - q(\bar{\psi})\mathbf{v}_d \cdot \nabla \vartheta + \frac{dq}{d\psi}(\bar{\psi})\hat{\psi} \frac{d\vartheta}{dt}_{\alpha_2}, \quad \chi(\alpha_2) = \int_0^t dt' [\dot{\alpha} - \langle \dot{\alpha} \rangle_{\alpha_2}] \quad (39) $$

Since $\delta_p \hat{\vartheta}(\alpha_2) + \chi(\alpha_2)$ is a periodic function of $\alpha_2$, $\omega_3 = \omega_d + \delta_p q(\bar{\psi})\omega_b$ and $\omega_d$ is well defined. This expression for $\omega_d$ is identical to the one used in [25], eq. 10. However, it has been chosen here to work with $\bar{\psi}$ rather than $\psi$, since the resonance is solely a function of the invariants, which simplify the analytical/numerical computation of $\delta W_K$.

Though, it should be noted that such a choice for $\bar{\psi}$ does not allow all the physics related the inertial layer to be retained, since the MHD displacement $\xi$ is naturally a function of $\psi$, not $\bar{\psi}$. Therefore, using a step function for $\xi$, inertial enhancement and global stability effects recovered in [25], where $\psi$ was chosen as radial variable, are not included in our model. However, for fast particles at low energies (~100keV), orbit width are relatively small, which means that considering $\xi(r) = \xi(\bar{r})$ should not impact significantly the solution of the kinetic eigenvalue problem.

The change of variable has been chosen as $P_\phi \rightarrow \bar{\psi} = -P_\phi/Ze$, $\hat{\psi} = mRv_\|/Ze$, such as $\bar{\psi}$ stands for the flux surface at the tips of the banana orbits for trapped particles. One should note that such a change of variable is arbitrary, allowable for all classes of particles and simplify expressions. The toroidal drift frequency therefore yields, at leading order in $\epsilon$

$$ \omega_d = \frac{q(\bar{r})E\lambda}{Ze\bar{r}R_0B_0} I_d, \quad I_d = \langle \cos \vartheta \rangle_{\alpha_2} + s(\bar{r}) \frac{m}{\epsilon\lambda E} \langle v_\| \rangle_{\alpha_2} \quad (40) $$

where $I_d$ is given in annex.

### 2.4. Derivation of $\delta W_K$

Using the energy principle, the kinetic contribution to (1) is given by

$$ \delta W_K = \frac{1}{2} \int_V d^3\mathbf{x} \xi^* \cdot \nabla \cdot \Pi_h \quad (41) $$
where $\Pi$ stands for the stress tensor. Considering the equation of motion for fast particles, discarding their inertia on the basis of their dilution in the whole plasma, $J_h \times B = \nabla \cdot \Pi_h$. Since $\xi = (b/B) \times \nabla \chi$, $\delta W_K$ can be rewritten as, using vectorial identities

$$ \delta W_K = -\frac{1}{2} \int d^3\mathbf{x} \mathbf{J}_{\perp h} \cdot \nabla \mathbf{X}^*, \quad \mathbf{J}_{\perp h} \equiv \int d^3\mathbf{v} F_{hes} \mathbf{v}_d \quad (42) $$

$J_{\perp h}$ being the perpendicular hot current. It should be noted here that the parallel component of the hot particles current has also to be taken into account, this term being incorporated into the functional $\delta W_{MHD}$ so that the current considered in Maxwell equations stands for the total current, as it should. Using the same normalization as in the previous part

$$ \lambda_K = \frac{1}{\xi_0 s_0 R_0 B_p} \int d^3 \mathbf{x} d^3 \mathbf{v} F_{hes} [e(\mathbf{v}_d \cdot \nabla)\mathbf{X}^*] \quad (43) $$

Therefore, using this expression, the resonant perturbation in Fourier space and applying Parseval's theorem

$$ \lambda_K = \frac{1}{s R B^2} \int_0^\infty \frac{\partial F_{eq}}{\partial E} \frac{\omega - \omega_*}{\sigma^2 E^2} \frac{\sigma^2 E^2}{R^2} \frac{1}{l_q} \frac{\omega + \delta \sigma (1-q)\omega - \omega}{\partial E} \quad (44) $$

The strategy that has been adopted to obtain precise numerical results for $\lambda_K$ is to apply analytically the strategy described in the main text. In the case of the resonant terms, the results are very准确 and the results depend on all integration variables. Therefore, computing the resonant integral analytically enhances the overall precision of $\lambda_K$. This strategy here differs from others linear codes such as [16] that use a grid adapted to the resonance patterns expected in the plane $(\psi, E)$. Therefore, in order to carry out such an analytical integration, one needs to express the resonant term as an unitary polynomial of $E$ or $v$.

$$ \lambda_K = \int_0^\infty \frac{\omega_*}{s_0 R_0 B_p} \frac{\sigma E^2 l_q^2}{\partial E} \frac{\partial F}{\partial R} \frac{\partial F}{\partial \hat{\omega}} \frac{1}{\partial \hat{r}} \frac{\hat{E} + \delta \sigma (1-q) \hat{\omega} \hat{E}^{1/2} - \hat{\omega}}{\partial \hat{r} \partial E} \quad (45) $$

where $\hat{E} = E/E_\alpha$, $\hat{v} = v/v_\alpha$, $\hat{\omega} = \omega/\tilde{\omega}_d$, $\tilde{\omega}_d = \omega_d[\tilde{r} = r_0, \lambda = (1 - \epsilon_0)^{-1}, E = E_\alpha] = E_\alpha/e_s B_0 r_0 R_0$, $r_0 = r_{q=1}$, $x = \tilde{r}/r_0$, $E_0 = x/q\sigma I_d$, $\tilde{\omega}_b = 2e_s B_0 r_0 / qm_\alpha v_\alpha I_{b,p}$, and where $\sigma_\parallel = \pm 1$ regarding co-passing/counter-passing particles. The interchange term reads

$$ \lambda_K^t = - \frac{1}{s_0 R_0 B_p} \int d^3 \mathbf{x} d^3 \mathbf{v} \frac{\partial F_{eq}}{\partial \tilde{r}} \frac{\sigma E}{R} e^{-i\theta} \quad (46) $$

One now needs an expression for the volume element in phase space using the coordinates $(\tilde{r}, \tilde{\vartheta}, \varphi, E, \lambda, \Phi)$, $\Phi$ being the gyro-angle. In toroidal geometry, assuming large aspect ratio and using $(v_\parallel, v_\perp, \Phi)$ coordinates

$$ \int d^3 \mathbf{x} = 2\pi R_0 \int_0^\infty \int_2\pi d\tilde{r} d\tilde{\vartheta}, \quad \int d^3 \mathbf{v} = 2\pi v_\perp dv_\parallel dv_\perp = 2\pi \int \sum_{E=1/2} \int_2^{1/2} \int_\lambda^{1/2} J dE d\lambda \quad (47) $$

where $\sigma_\parallel = \text{sgn}(v_\parallel)$, $J$ is the Jacobian of the change of variable $(v_\parallel, v_\perp) \to (E, \lambda)$. Since $v_\parallel = \sigma_\parallel(2E/m)^{1/2}(1 - \lambda/H)^{1/2}$, $v_\perp = (2\lambda E/mH)^{1/2}$, one has the following derivatives

$$ \frac{\partial v_\parallel}{\partial E} = \sigma_\parallel (2Em)^{-1/2} (1 - \lambda/H)^{1/2}, \quad \frac{\partial v_\parallel}{\partial \lambda} = -\sigma_\parallel (2Em)^{1/2} (2H)^{-1} (1 - \lambda/H)^{-1/2} \quad (48) $$

$$ \frac{\partial v_\perp}{\partial E} = (\lambda/2mEH)^{1/2}, \quad \frac{\partial v_\perp}{\partial \lambda} = (E/2m\lambda H)^{1/2} \quad (49) $$
Therefore, using the following identity, phase space in our set of variables is expressed as

$$ \int_{-\pi}^{\pi} \int_{0}^{H} \int_{0}^{(1-\epsilon)^{-1}} \int_{-\vartheta_0}^{\vartheta_0} d\vartheta \quad d\lambda = \quad d\lambda \quad d\vartheta \quad (51) $$

$$ \int d^3\mathbf{x} d^3\mathbf{v} = \sum_{\sigma_{\parallel}=\pm 1} \pi^2 R_0 \frac{2}{\sqrt{3}} \int_{0}^{3/2} \int_{\bar{r}d\bar{r}} r_0 \int_{(1-\epsilon)^{-1}}^{\epsilon} d\lambda \int_{-\vartheta_0}^{\vartheta_0} 1 - \frac{\lambda}{\bar{\lambda}}^{-1/2} d\vartheta \int_{\vartheta_0}^{\infty} E^{1/2} dE \quad (52) $$

where the sum over the parallel velocity is two for trapped particles. Regarding the different terms of $\lambda_K$, considering only even terms in $\vartheta$ and the bounce-averaging formalism developed, one gets

$$ \lambda_K^{\text{res}} = \frac{2\pi^3\epsilon}{s_0 r_0 B_0} \frac{2}{m} \sum_{p,0} \int_0^{3/2} \int_{\bar{r}d\bar{r}} (1-\epsilon)^{-1} d\lambda \sigma_{\parallel} I_b \frac{l_q^2}{l} \int_0^{\infty} E^{5/2} \frac{\partial \bar{r} F_{\text{eq}}/E\alpha - \bar{r} \hat{\omega} \partial_E F_{\text{eq}}/qr}{\hat{E} + \delta \sigma (1-q)E |\hat{\Omega}|} dE \quad (53) $$

$$ \lambda_K^{\text{int}} = - \frac{4\pi^3\epsilon_0}{s_0 r_0 B_{p,0}} \frac{2}{m} \sum_{\sigma_{\parallel}=\pm 1} \int_0^{3/2} \int_{\bar{r}d\bar{r}} (1-\epsilon)^{-1} d\lambda \sigma_{\parallel} I_b \int_0^{p_{\parallel}} E^{3/2} \frac{\partial F_{\text{eq}}}{\partial \bar{r}} dE \quad (54) $$

where $l_c = \langle \cos \vartheta \rangle_{\alpha_2}$, $l_b = \langle 1 \rangle_{\alpha_2}$. After little algebra, one can show that these two expressions are identical to $\delta W_1$ and $\delta W_2$ [14]. This re-derivation has been done to obtain expressions suited for analytical/numerical computation. Since an analytical treatment is preferred for the energy integral, calculations have to be specific to a given equilibrium particles distribution function, which is done for an isotropic slowing-down distribution function in Appendix B. Now that the kinetic term has been derived, let us discuss the resolution of the dispersion relation

### 3. Resolution of the dispersion relation with MHD-K and linear results

#### 3.1. MHD-K

The linear code MHD-K solves equation (1) using a partial analytical approach. $\lambda_K$ is integrated over all phase space, however, as it was shown, its integrand only depends on the triplet $(\bar{r}, E, \lambda)$. As shown in Appendix B, given an equilibrium function for fast particles, one can solve analytically the energy integral. Such an analytical treatment of the resonance brings a great deal of precision in the computation of $\lambda_K$. This is due to the fact that a 3D grid would have been inaccurate, the poles of the integrand being multiple and dependent on the triplet $(\bar{r}, E, \lambda)$. The other integrals over the pitch-angle and radial position are solved numerically using elliptic integrals of the first and second kind, as presented in Appendix A.

The fluid term $\lambda_H = \gamma_{MHD} \tau_A$ of the dispersion relation is calculated using the coupled codes CHEASE [26], that computes the 3D toroidal MHD equilibrium and XTOR-2F [27], that time iterates the equilibrium through the bi-fluid resistive extended MHD equations. All equilibrium quantities and profiles used in MHD-K are directly taken from CHEASE. It is important to note that $\gamma_{MHD}$ is the linear growth rate computed by XTOR-2F for a given equilibrium in absence of fast particles. Therefore, it is assumed that the fast particle current is negligible regarding the total current. Fluid contributions are therefore not counted twice, since the total current response is present in $\delta W_{MHD}$, and the fast particles' pressure response in $\delta W_K$ [14].

MHD-K assumes, as does the linear model, that the equilibrium has concentric circular magnetic flux surfaces without Shafranov shift. Therefore, CHEASE is used with circular plasma cross section. However with circular flux surfaces, a Shafranov shift is always present when solving Grad-Shafranov equation. This is handled in MHD-K by shifting its geometric axis to the magnetic axis of CHEASE, assuming a larger major radius than expected.

Then, (1) is solved through a method somewhat simpler than in [28]. Since the goal of the linear code is to verify non-linear hybrid kinetic-MHD codes, it is not necessary to use a complex algorithm locating the zeros of an analytic function. Instead, the value of the kinetic term $\lambda_K$ is computed on a 2D grid
in the complex plane $(\omega, \gamma)$, then $1/|D(n_{h,0}, \omega + i\gamma)|^2$ is computed all over the grid, its maximum being considered as the solution of (1) for a given kinetic density on the magnetic axis $n_{h,0}$.
By taking the maximum of $1/|D(n_{h,0}, \omega + i\gamma)|^2$ as the solution of (1), it is assumed that the solutions are unique, which is not necessarily the case. In fact, it has been observed that for the internal kink in presence of fast particles, two branches co-exist. The kink branch at low fast particles density, and the fishbone branch at high density. But since the kink and fishbone branches are stabilized respectively at high and low density, only one branch exists in the unstable region $\gamma > 0$. Only the unstable solutions are of interest, since the goal of MHD-K is to be used as a tool to verify hybrid codes, where only unstable modes can exist.

## 3.2. Effects of passing particles on the fishbone mode

In the following, linear results obtained from MHD-K are presented. These results have been obtained from an ITER-like equilibrium with $R_0 = 6.2m, B_0 = 5.3T, a = 1.86m, S = 1.10^7, \omega_{*i} = 0$, parabolic-like q profile with $q_0 = 0.8$ on the magnetic axis. Bulk and fast particles profiles have the form $P = P_0(1 - r^\alpha)^\beta$, where for the bulk density $n_i, n_{i,0} = 2.10^{19}m^{-3}, \alpha = \beta = 2$, for the bulk temperature $T_i, T_{i,0} = 30keV, \alpha = \beta = 2$ and for the fast particles density $n_h, n_{h,0}$ is a parameter and $\alpha = 2, \beta = 6$.

The fast particles are distributed along an isotropic slowing-down equilibrium distribution function, as presented in Appendix B. The peak energy of the particle is set to $E_\alpha = 3.5MeV$, which is the energy of fusion-born alphas. In later comparison with hybrid codes, this peak energy might be decreased for precise comparisons between the linear model and numerical results. At high energy $\gamma \approx 1MeV$, particles follow non-standard orbits known as potato orbits [29], as the particle's orbit width increases with its energy. However, it has been assumed in the linear model that particles follows thin banana orbits, therefore lower energies $\sim 100keV$ would be preferred.

By varying $n_{h,0}$ from 0 to $2.10^{17}m^{-3}$, results shown in Figure 1 have been obtained. Considering

<figure><img src="image_3.png" /><figcaption>Figure 1: Linear growth rate of the kinetic kink instability, $\beta_\alpha$ the fast particles beta on the magnetic axis, considering respectively all particles, only passing particles and only trapped particles</figcaption></figure>

the full model curve in Figure 1, one can see that the typical behavior for the kinetic kink instability is retrieved. At low fast particles density, the internal kink is stabilized and then, past a threshold in $\beta_\alpha$, $\beta_\alpha/\beta \sim 0.35$ here, the fishbone mode is dominant and the kink instability is destabilized.

The variation of $\gamma$ regarding $n_{h0}$ needs to be explained here. This variation is due to two competing physical effects. By looking closely at (1), assuming $\omega_{*i} = 0$, $\gamma\tau_A = \gamma_{MHD}\tau_A + n_{h0}\text{Re}[\tilde{\lambda}_K(\omega + i\gamma)]$, $\tilde{\lambda}_K = \lambda_K/n_{h0}$. Only a variation of the real part of $\lambda_K$ can explain the variation of $\gamma$. As it can be seen on Figure 2, the overall (trapped and passing particles) fluid term inside $\text{Re}[\lambda_K]$ is due to the negative density gradient and the mostly positive sign of $l_c$ in the $(\psi, \lambda)$ plane. Therefore it is not the one producing the stabilization. The real resonant term $\text{Re}[\lambda_{\text{res}}]$ is mostly negative in phase space
according to the sign of $I_d$. Then, the stabilization is produced by the non-reactive part, a fluid-like effect of $Re[\lambda^{res}]$, since at low frequency it can be shown that in most of the phase space $Re(\omega) = n \Omega$. When the fast particles density is increased, this fluid-like effect competes with the resonant effect, because there are now significant zones of the phase space where $Re(\omega) = n \Omega$, leading to the destabilization. This variation of $Re[\lambda^{res}]$ is illustrated in Figure 2. Physically, it can be explained by the fact that at low density, the fast particles behave mostly as a fluid since few particles occupy the small regions in which they interact resonantly with the MHD mode. As the density is increased, more particles occupy the resonant zones, leading to a prevalence of the resonant effect on the fluid effect.

It can also be seen on Figure 1 that passing particles have a significant effect on the fishbone mode, since the results of the full model are not retrieved when only considering trapped particles. This can firstly be explained by the fact that fluid terms from passing particles are of the same order of magnitude as for the trapped particles, as shown in Figure 2. Secondly, as it will be shown in the following part, the MHD mode also reacts resonantly with passing particles. This is also illustrated

<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of fluid and real resonant part of the kinetic term, underlying the competitive physical effects</figcaption></figure>

by Figure 3, where the pulsation of the MHD mode is presented for the three models. When solely considering passing or trapped particles, the MHD mode frequency is not as high as it is when all classes of particles are considered. This due to the fact that the kink instability has less resonant interactions with the energetic particles, which implies that resonant surfaces also exist inside the passing part of phase space.

### 3.3. Resonances plane in phase space

In Figure 4, the resonance curves at a fixed radial position $\bar{r} = 0.8 r_{q=1}$ are displayed. They are obtained by plotting the integrand of $\delta W_K$ in the 3D space $(\bar{r}, E, \lambda)$. Regarding (20), only two particle modes contribute here, $n = (0, 1, 1)$ for passing particles and $n = (0, 0, 1)$ for trapped particles. In the phase space zone concerned by passing particles, two different resonant curves, associated with the same particle mode, can be seen. They correspond to the co-passing and counter passing contributions to the mode. Regarding the amplitude of the curves, in most cases, the trapped curve has an amplitude up to five times higher than passing curves. However, passing curves usually cover a broader area in phase space. Therefore, this result underlines the necessity to consider passing particles in the linear model to properly describe the kinetic internal kink mode, since depending on the imposed equilibrium function in phase space, the resonant effects could be due entirely to passing particles.
<figure><img src="image_5.png" /><figcaption>Figure 3: Linear pulsation of the kinetic kink instability, $\beta_\alpha$ the fast particles beta on the magnetix axis, considering respectively all particles, only passing particles and only trapped particles</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Resonant curves in the $(E, \lambda)$ plane, for a fixed radial position $\bar{r} = 0.8r_0$, $r_0$ being the $q = 1$ radius</figcaption></figure>

**4. Conclusion**

A new linear model, including more physical effects than some previous ones, has been derived here using a Hamiltonian angle-action formalism. On its basis, the linear code MHD-K has been developed, solving non-perturbatively the kinetic internal kink dispersion relation. The linear results obtained from MHD-K have shown the relevance of considering passing particles in linear models, since they are also a resonant drive to the emergence of the fishbone mode. This new linear tool aims at testing non-linear hybrid kinetic codes during their linear phase, in order to verify their correct implementation. Such a verification is currently underway with the code XTOR-K [7][8], a hybrid kinetic-MHD using for the MHD part the solver from [27], and for the kinetic part a 6D PIC module.

This work has been carried out within the framework of the French Research Federation for Fusion Studies, and of the EUROfusion Consortium. It has received funding from the Euratom research
and training programme 2014-2018 under grant agreement No 633053 for the project WP17-ENR-CEA-06. This work is part of the project AMICI funded by the Agence Nationale pour la Recherche (ANR-14-CE32-0004-01). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## Appendix A. Bounce-averaging formalism

In this section, bounce-averaging calculations are detailed for both passing and trapped particles. Any bounce-averaged quantity $<F>$ is calculated through

$$ <F> = \int_{C v_{\parallel}}^{\int_{C v_{\parallel}}^{F} dl / \int_{C v_{\parallel}}^{dl}} \frac{2E^{1/2}}{m} \left(1 - \frac{\lambda}{H}\right)^{1/2} \quad (A.1) $$

where C is the path followed by the particle. The parallel velocity can be expressed as follows, considering $E = 1/2mv^2$, $\mu = mv^2/2B$, $\lambda = \mu B_0/E$, $H = R/R_0$

$$ v_{\parallel} = \frac{2E^{1/2}}{m} \left(1 - \frac{\lambda}{H}\right)^{1/2} \quad (A.2) $$

Assuming circular flux surfaces, it yields, for $\vartheta_0$ the bounce angle for trapped particles and $\vartheta_0 = \pi$ for passing ones

$$ <F> = \int_{-\vartheta_0}^{\vartheta_0} F \left(1 - \frac{\lambda}{H}\right)^{-1/2} \frac{d\vartheta}{2\pi} / \int_{-\vartheta_0}^{\vartheta_0} \left(1 - \frac{\lambda}{H}\right)^{-1/2} \frac{d\vartheta}{2\pi} \quad (A.3) $$

using $\epsilon = r/R_0$, one has

$$ \frac{\lambda}{H} = 1 - \frac{1}{1 - \lambda + \epsilon\lambda} \quad 1 - \frac{\sin^2(\vartheta/2)}{y^2} \quad (A.4) $$

with $y^2$ being expressed as

$$ y^2(x, \lambda) = \frac{1 - \lambda + \epsilon\lambda}{2\epsilon\lambda} \quad (A.5) $$

It can be noted that [12] uses a different definition for $y^2$. One can distinguish passing from trapped particles by a simple condition on $y^2$. Trapped particles are characterized by the existence of two turning points in the poloidal plane where their parallel velocity is null. Therefore, the condition on having trapped particles translates itself, according to (A.2), as $(1 - \lambda/H)^{1/2} = 0$, which implies, assuming $\lambda < 1/1 - \epsilon$, that $0 < y^2 < 1$. Passing particles are therefore characterized by $1 < y^2 < +\infty$.

From these conditions on $y^2$, one can deduce the corresponding conditions on the pitch angle $\lambda$. Respectively, for $y^2 = 0, 1, +\infty$, the corresponding $\lambda$ are $\lambda = 1/1 - \epsilon, \lambda = 1/1 + \epsilon, \lambda = 0$. Therefore, for trapped particles, $1/1 + \epsilon < \lambda < 1/1 - \epsilon$, and for passing particles, $0 < \lambda < 1/1 + \epsilon$.

These expressions enable to write $<F>$ as

$$ <F> = \int_{-\vartheta_0}^{\vartheta_0} F \left(1 - y^{-2} \sin^2(\vartheta/2)\right)^{-1/2} \frac{d\vartheta}{2\pi} / \int_{-\vartheta_0}^{\vartheta_0} \left(1 - y^{-2} \sin^2(\vartheta/2)\right)^{-1/2} \frac{d\vartheta}{2\pi} \quad (A.6) $$

Let us now derive expressions for $<\cos\vartheta>$, $<\vartheta \sin\vartheta>$, $<v_{\parallel}^2>$ and $<\cos(q\vartheta)>$ for both passing and trapped particles, which are appearing in the derivation of $\lambda_k$.

### Appendix A.1. Passing particles

For passing particles, $\vartheta_0 = \pi$, the frequency to perform a poloidal revolution is defined as

$$ \frac{\omega_b}{2\pi} = \frac{(2E/m)^{1/2}}{qR_0 I_{b,p}} \quad (A.7) $$
Moreover, knowing that $\cos \vartheta = 2(1 - y^2 \sin^2 u) - 1$

$$ \int_{-\vartheta_0}^{\vartheta_0} \cos \vartheta \left(1 - \frac{\lambda}{H}^{-1/2} d\vartheta = 4 \right) 2\lambda \frac{r}{R_0} \frac{-1/2}{0} \int_{0}^{\pi/2} \frac{\cos \vartheta}{\cos(\vartheta/2)} d\vartheta = 4 \right) 2\lambda \frac{r}{R_0} \frac{-1/2}{0} \left[ 2E(y^2) - K(y^2) \right] \quad (A.19) $$

which gives

$$ \langle \cos \vartheta \rangle_{\tau} = I_{c,t} = 2 \frac{E(y^2)}{K(y^2)} \quad (A.20) $$

Performing the same integration by parts as before for $\langle \vartheta \sin \vartheta \rangle_{\tau}$, where this time the first term is dropped since $1 - y^{-2} \sin^2 \vartheta_0/2 = 0$

$$ \int_{-\vartheta_0}^{\vartheta_0} \vartheta \sin \vartheta \left(1 - \frac{\lambda}{H}^{-1/2} d\vartheta = 4 \right) 2\lambda \frac{r}{R_0} \frac{-1/2}{4y} \int_{-\vartheta_0}^{\vartheta_0} \left(1 - y^{-2} \sin^2(\vartheta/2)\right)^{1/2} d\vartheta \quad (A.21) $$

Therefore, noticing that $\cos^2 u = y^{-2}(1 - y^2 \sin^2 u) + 1 - y^{-2}$

$$ \langle \vartheta \sin \vartheta \rangle_{\tau} = I_{s,t} = \frac{4y^2}{K(y^2)} \int_0^{\pi/2} \frac{\cos^2 u}{(1 - y^2 \sin^2 u)^{1/2}} du = 4 \frac{E(y^2)}{K(y^2)} + 4(y - 1) \quad (A.22) $$

For $\langle v^2 \rangle$, it can be shown easily, using the same expression for $\cos^2 u$, that

$$ \langle v^2 \rangle_{\parallel} = \frac{E\lambda\epsilon}{m} \langle \vartheta \sin \vartheta \rangle \quad (A.23) $$

The expression of $\langle \cos q\vartheta \rangle$ is defined identically for both passing and trapped particles.

## Appendix B. Kinetic contribution using a slowing down distribution

Using a normalized slowing-down distribution function, different from [12]

$$ F_{eq}(\bar{r}, v) = \frac{3 m^{3/2}}{4\pi \ln[1 + (v_{\alpha}/v_c)^3]} \quad (B.1) $$

with $v_c$ the critical speed at which fast particles give out as much energy to bulk ions and electrons, expressed as, with $m_h$ the mass of hot particles

$$ v_c^3(\bar{r}) = \frac{3 \sqrt{\frac{\pi}{m_e}} m_e}{4 m_h} \frac{2T_e(\bar{r})}{m} \quad (B.2) $$

Considering the following change of variable $\hat{E} = \hat{v}^2$ in order to express the resonance as a second order polynomial, it yields, for the main contribution to $\lambda^{\text{res}}$ proportional to $\omega_*$, using $v_1 = (1 - q)|\hat{\Omega}_2|$

$$ \lambda^{\text{res},\omega_*} = \frac{3\pi^2 \epsilon_0 E_\alpha}{p,0} \sum_{\kappa} \int_{r_0}^{r_0} \frac{dn_\alpha}{d\bar{r}} \int_{(1-\epsilon)^{-1}}^{1} d\lambda \sigma_l \int_{q}^{l^2} \frac{1}{\hat{v}^3} d\hat{v} \\ \qquad \kappa \quad 2s_0 r_0 B^2 \ln[1+(v_\alpha/v_c)^3] \sigma_{\parallel=\pm 1} \int_0^0 d\bar{r} \quad (B.3) $$

It should be noted that from the term $1/(v^3 + \zeta^3)$, only the contribution $1/v^3$ has been kept. This is a reasonable approximation for highly energetic particles, such as fusion born alphas, since in this case $(v_c/v_\alpha)^3 \sim 10^{-1}$. For fast particles with lower energies ($\sim 100 \text{keV}$), the reverse needs to be done.

One now needs to identify the resonant roots, that are functions of the considered $(\bar{r}, \lambda)$ couple. Given the usual determinant $\Delta = E_0[\delta_p E_0 v_1^2 + 4\hat{v}]$, the roots yields

$$ \hat{v}_{\pm}(r, \lambda) = \frac{-\delta_p \sigma_{\parallel} v_1(\bar{r}, \lambda) E_0(\bar{r}, \lambda) \pm \Delta(\bar{r}, \lambda)}{2} $$
Valon the general case, the roots are complex, since $|v_1E_0| \sim 1$ for $1-q \sim 10^{-1}$, with opposed IOP Publishing
IOP parts should be noted that for CTE R likes, equilibrio and 0.03-0.9 doi:10.1088/17426596/4/2/1082003001
particular case of trapped particles, the roots are simply given by $\pm \sqrt{E_0\bar{\omega}}$. Then, using these roots, one can perform a decomposition in simple elements, which yields, after straightforward algebra

$$ \int_0^1 \frac{2\hat{v}^3}{\hat{v}^2 + \delta \sigma E v \hat{v} - E \hat{\omega}} d\hat{v} = \left(\hat{v} - \hat{v}_-\right) \hat{v}_- + \left(\hat{v}_- - \hat{v} + \delta \sigma E v \hat{v} - E \hat{\omega}\right) d\hat{v} \quad (B.5) $$

which expands to

$$ \int_0^1 \frac{\hat{v}^2 + \delta \sigma E v \hat{v} - E \hat{\omega}}{p \parallel 0 1 0} d\hat{v} = \left(\hat{v} - \hat{v}_-\right) \times \hat{v}_- \frac{1}{2} + \hat{v}_- \left(1 + \hat{v}_- \ln 1 - \frac{1}{\hat{v}}\right) \hat{v}_- + \left(\hat{v}_- - \hat{v} + \frac{1}{2} + \hat{v}_+ + \hat{v}_+ \ln 1 - \frac{1}{\hat{v}}\right) \hat{v}_+ \quad (B.6) $$

When computing this integral, one needs to handle carefully the complex logarithm cut-off. When varying $(\bar{r}, \lambda)$ to solve the remaining integrals, the resonant roots need to vary continuously. For a cut-off located on the negative part of the real axis in the complex plane, it can be shown that the cut off happens when $Im(v_\pm)$ switches sign while $Re(1 - 1/v_\pm) < 0$. Since the complex logarithm is given by $\ln_c z = \ln|z| + 2i \arctan(Im(z)/[Re(z) + i\bar{\omega}])$, a factor $2i\pi$ is added to the complex logarithm when the poles are crossing the cut off half axis, the sign depending on the crossing direction.

The resonant term has been divided in two parts regarding the diamagnetic frequency and the mode frequency. The contribution proportional to the mode frequency is supposed to be lower than the other, but it is derived here anyway

$$ \lambda_K^{\text{res}, \hat{\omega}} = \sum_{\sigma \parallel \pm 1} \int_{r_0}^{r_0} \frac{\chi^2 n}{q^\alpha \bar{r}^\alpha} \int_{0}^{(1-\epsilon)^{-1}} \frac{l^2}{d\lambda \sigma l_b^q l^q} \\ \times \frac{1}{1 + \delta \sigma E v - E \bar{\omega}} + \frac{3}{2} \int_{0}^{1} \frac{\hat{v}}{\hat{v}^2 + \delta \sigma E v - \hat{v} - E \bar{\omega}} d\hat{v} \quad (B.7) $$

with

$$ \int_0^1 \frac{\hat{v} d\hat{v}}{p \parallel 0 1 0} = \frac{1}{\hat{v}^2 + \delta \sigma E v \hat{v} - E \hat{\omega}} = \hat{v} - \hat{v}_- \times \hat{v}_- \ln(1 - 1/\hat{v}_-) - \hat{v}_+ \ln(1 - 1/\hat{v}_+) \quad (B.8) $$

The interchange term simply reads

$$ \lambda_K^{\text{int}} = - \frac{3\pi^2 \epsilon_0 E_\alpha}{s_0 r_0 B_{p0} \ln[1 + (v_\alpha/v_c)^3]} \int_0^1 \frac{dn}{d\bar{r}} \int_{0}^{(1-\epsilon)^{-1}} \frac{l^2}{l^2} \times \frac{\alpha \bar{r} d\bar{r}}{\sigma l_b^q d\lambda} \quad (B.9) $$

[1] K. McGuire and al, PRL, 50, 12 (1983)

[2] M. Nave and al, Nuclear Fusion, Vol.31, No.4 (1991)

[3] M. J. Mantsinen and al, Plasma Phys. Control. Fusion 42 (2000) 1291-1308

[4] Y. Pei, N. Xiang, Y. Hu, Y. Todo, G. Li, W. Shen, L. Xu, Physics of Plasmas 24, 032507 (2017)

[5] G. Y. Fu and al, Physics of Plasmas 13, 052517 (2006)

[6] F. Wang, G. Y. Fu, J. A. Breslau, J. Y. Liu, Physics of Plasmas 20, 102506 (2013)

[7] D. Leblond, PhD thesis, (2011)

[8] H. Lutjens and al, Description of the hybrid code XTOR-K, to be submitted (2018)

[9] D. Edery, X. Garbet, J-P. Roubin, A. Samain (1992), PPCF, 34, 6, 1089-1112

[10] L. Chen, R. B. White, M. N. Rosenbluth, PRL, 52, 13, (1984)

[11] R. B. White, F. Romanelli and M. N. Bussac, Phys. Fluids B 2, 745 (1990); doi: 10.1063/1.859312

[12] Coppi, B., S. Migliuolo, F. Pegoraro and F. Porcelli (1990). Phys. Fluids B2, 927

[13] F. Porcelli, PPCF, 33, 13 1601-1620, (1991)

[14] F. Porcelli, R. Stankiewicz, W. Kerner, H.L. Berk (1994) Physics of Plasmas, 1,470

[15] C. Z. Cheng Physics Reports 211, No. 1 (1992) 1-51. North-Holland
[16] F. Nabais and al. Plasma Science and Technology, Vol.17, No.2, Feb. 2015

[17] Freidberg, J. (2014). Ideal MHD. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511795046
Vainio, P., Copi, R. J., Hastie, S., Migliuolo, F., Pegoraro, and F. Porcelli. Phys. Lett. A 32, 267 (1988).

[19] F. Zonca, L. Chen, Physics of Plasmas 21, 072120 (2014)
[20] F. Zonca, L. Chen, Physics of Plasmas 21, 072121 (2014)

[21] X. Garbet, PhD thesis, Annex A, 1988

[22] A.N. Kaufman, Phys. FLuids, 15, 1063 (1972)

[23] H.E. Mynick and J.A. Krommes, Phys. Rev. Lett. 43, 1506 (1979)

[24] A.J.Brizard, T.S.Hahm, Reviews of modern physics, 79, 2007

[25] J. Graves, Plasma Phys. Control. Fusion 55 (2013) 074009

[26] H. Lütjens, A. Bondeson, O. Sauter, Computer Physics Communications 97 (1996) 219-260

[27] H. Lütjens, J-F. Luciani, JCP, (2010), 229, 8130-8143

[28] B. Davies, Journal of Computational Physics 66, 3649 (1986)

[29] F. Porcelli, R. Stankiewicz, H. L. Berk, Y. Z. Zhang, Phys. FLuid B 4 (10), 1992