# Global linear stability analysis of kinetic Trapped Ion Mode (TIM) in tokamak plasma using spectral method

Debraj Mandal, M. Lesur, E Gravier, J.-N. Sama, A. Guillevic, Y. Sarazin, X. Garbet

► To cite this version:

Debraj Mandal, M. Lesur, E Gravier, J.-N. Sama, A. Guillevic, et al., Global linear stability analysis of kinetic Trapped Ion Mode (TIM) in tokamak plasma using spectral method. Plasma Physics and Controlled Fusion, 2023, 65 (5), pp.055001. 10.1088/1361-6587/acbe64. hal-03864496v2

HAL Id: hal-03864496

https://hal.science/hal-03864496v2

Submitted on 15 Mar 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.

Public Domain
# Global linear stability analysis of kinetic Trapped Ion Mode (TIM) in tokamak plasma using spectral method

D. Mandal,¹ M. Lesur,¹ E. Gravier,¹ J. N. Sama,¹ A. Guillevic,¹ Y. Sarazin,² and X. Garbet²

¹Université de Lorraine, CNRS, IJL, F-54000 Nancy, France

²Association Euratom-CEA, CEA/DSM/DRFC, Centre de Cadarache, F-13108, St-Paul-Lez-Durance, France

Trapped ion modes (TIM) belong to the family of ion temperature gradient (ITG) modes, which are one of the important ingredients in heat turbulent transport at the ion scale in tokamak plasmas. A global linear analysis of a reduced gyro-bounce kinetic model for trapped particle modes is performed, and a spectral method is proposed to solve the dispersion relation. Importantly, the radial profile of the particle drift velocity is taken into account in the linear analysis by considering the magnetic flux $\psi$ dependency of the equilibrium Hamiltonian $H_{eq}(\psi)$ in the quasi-neutrality equation and equilibrium gyro-bounce averaged distribution function $F_{eq}$. Using this spectral method, linear growth-rates of TIM instability in presence of different temperature profiles and precession frequencies of trapped ions, with an approximated constant Hamiltonian and the exact $\psi$ dependent equilibrium Hamiltonian, are investigated. The growth-rate depends on the logarithmic gradient of temperature $\kappa_T$, density $\kappa_n$ and equilibrium Hamiltonian $\kappa_\Lambda$. With the exact $\psi$ dependent Hamiltonian, the growth rates and potential profiles are modified significantly, compared to the cases with approximated constant Hamiltonian. All the results from the global linear analysis agree with a semi-Lagrangian based linear Vlasov solver with a good accuracy. This spectral method is very fast and requires very less computation resources compared to a linear version of Vlasov-solver based on a semi-Lagrangian scheme.

## I. INTRODUCTION

Low frequency and low wavenumber turbulence, which is mainly generated by ion temperature gradient (ITG) [1–3] and trapped electron mode (TEM) [4, 5] instabilities, plays a dominant role in the anomalous radial energy and particle transport in magnetically confined fusion plasmas. ITG is an important ingredient in anomalous ion heat transport in tokamak whereas TEM turbulence drives electron particle and heat transport. Ion temperature gradient driven modes are frequently observed and relevant in tokamak plasma experiments [6, 7]. Trapped ion mode (TIM) belongs to this family of ion temperature gradient modes, and is driven by the resonant motion of trapped ions [8]. The trapped ion instability is characterized by frequencies of the order of the trapped ion precession frequency and radial scales of the order of several banana widths. Though TIM modes have lower frequency and long wavelength compared to TEM modes, TIM has similar physical mechanisms as TEM in the linear regime. Both TIM and TEM are driven by precession frequency of ions and electrons respectively, and are generated in presence of ion and electron equilibrium gradients respectively [9, 10]. However, in the nonlinear regime they depart from each other due to different response of zonal flows [11–13]. It is essential to properly estimate the linear growth rate of these trapped-particle instabilities to understand their influence on nonlinear saturation, the turbulent nature of the system, and the associated transport. Although TEM is more directly relevant to tokamak turbulence and transport, TIM in numerical simulation is much more tractable and shares essential similarities with TEM in linear phase. Therefore here we focus on TIM instability.

During last few decades, ion turbulence in magnetically confined plasmas has been intensively studied in fluid simulations [14–16], fluid-kinetic hybrid electron model for studying low frequency electro-magnetic turbulence [17], particle-in-cell based gyrokinetic simulations [18–20], continuum Vlasov approach gyrokinetic simulations in Eulerian grid [21], $\delta f$-particle-in-cell based simulations for bounce-averaged kinetic equations obtained by phase-space Lagrangian Lie-perturbation theory [22, 23], and semi-Lagrangian based reduced bounce averaged gyrokinetic simulations [8, 9, 24]. However, a kinetic model is necessary to capture the kinetic features of trapped ion modes, and a full 5-dimension (5D) gyrokinetic model demands more computer resources, we consider a reduced 4-dimension (4D) gyro-kinetic model with averaging over gyro-motion and the banana-orbit motion, and adiabatic response of passing particles [8]. This allows us to study the trapped ion mode (TIM) instability, driven through the resonant interactions of trapped ions with a wave, separately from the ITG instability which is generated in presence of ion pressure gradient.

The trapped ion mode instability is driven by the ion temperature gradient $\nabla_{\perp} T_i$ and $\nabla B$ drift. In the unfavorable curvature region of a tokamak, temperature gradient is aligned with the magnetic field gradient, particles in the lower temperature region drift more slowly than in the higher temperature region, which yields a charge separation in the presence of a perturbed density profile. The electric field generated from this charge separation creates a $E \times B$ drift motion, which enhances the initial perturbation, hence leading to an instability. Though the nonlinear evolution of this TIM instability is interesting in the context of turbulent transport, it is first worth studying the linear properties of this instability. A
local-linear analysis of TIM instability using this reduced gyro-kinetic model was reported by Drouot et. al. [9], where the values of all relevant parameters were considered at a particular radial location. The linear version of the simulation code TERESA [8, 25, 26], which is based on a semi-Lagrangian method, is able to solve the modeled equations in the linear limit with radial profiles of all relevant parameters [27]. Here we are presenting an alternative global-linear analysis by solving the modeled equations of this reduced gyro-kinetic model in the linear limit using a spectral method [28, 29], which takes into account the radial profiles of all relevant parameters, and saving a lot of computational resources compared to the linear version of TERESA code.

Moreover, till now, the trapped ion mode instability, and its growth-rate in the presence of different temperature profiles and variations in precession frequency have been studied with the assumptions that the temperature profile varies linearly with $\psi$ and that the precession frequency remains constant throughout the simulation box, which is a good approximation for a simulation region, situated sufficiently far from the edge (i.e. the last closed flux surface (LCFS) in our model) and the core region of a tokamak. Here using this solver, based on a spectral method, we have studied the effects on the TIM instability of different temperature profiles $T_i(\psi)$, and precession frequency profiles $\Omega_D(\psi, \kappa)$ as a function of the magnetic flux surface $\psi$ and trapping parameter $\kappa$, which is associated with the pitch angle variation of trapped particles. All the temperature and precession frequency profiles are relevant to different experimental observations. Depending on the temperature profiles and precession frequency values we found different growth rates of the TIM instability and different potential solutions of the system, which are significantly different from the previous local linear analysis results [9]. Additionally for incorporating the variation of particles drift velocity with magnetic flux $\psi$, we have modified the previous gyro-bounce averaged kinetic model by introducing the $\psi$ dependency of the equilibrium Hamiltonian $H_{eq}$ in the quasi-neutrality equation and the equilibrium distribution function, which was neglected in the previous model where the equilibrium Hamiltonian was approximated by the energy $E$ of the particles: $H_{eq} \approx E$ [25]. This new modification changes the growth-rate profile of TIM mode instability, and changes the potential solution of the system. The growth-rate depends on the logarithmic gradients of temperature, density and equilibrium Hamiltonian, respectively $\kappa_T$, $\kappa_n$ and $\kappa_\Lambda$. All the results from this spectral method are successfully compared with the new linear version of TERESA, which incorporates the $\psi$ and $\kappa$ dependency of precession frequency. Therefore this newly proposed method for global linear analysis of TIM instability, and the associated results help to understand the threshold/beginning of turbulent transport in tokamak by TIM modes, and also explain possible dominant mode of this instability in nonlinear regime.

Our paper is organized as follows, sec. II presents the

bounce-averaged gyrokinetic model, where a brief description of the previous nonlinear model, along with crucial upgrades are discussed. Also the formulation of global linear analysis of this nonlinear model is presented in this section. Sec. III presents a numerical solver based on spectral method for solving the dispersion relation, which is derived from global linear analysis. The effect of different temperature profiles and precession frequency profiles on the TIM instability in the limit $H_{eq} \approx E$ is discussed in sec. IV. In sec. V, the effect of the inverse gradient length of equilibrium Hamiltonian $\kappa_\Lambda$, on the TIM instability is presented. Sec. VI presents the conclusions.

## II. THE BOUNCE-AVERAGED GYROKINETIC MODEL AND ITS MODIFICATIONS DUE TO EXACT HAMILTONIAN ($H(\psi, \kappa)$)

An electrostatic reduced collisionless bounce-averaged gyrokinetic model was developed by Depret, Sarazin and Darmet [8, 9, 25, 26]. We adopt that model to study the stability of trapped ion modes. In this model the system evolves on a timescale of the order of the trapped particle precession frequency $\omega_D$, which allows to filter out the large frequencies $\omega_c$ (cyclotron frequency) and $\omega_b$ (bounce frequency around banana orbit) ($\omega_D \ll \omega_b \ll \omega_c$), and to simplify the effect of length scales $\rho_c$ (gyro-radius) and $\delta_b$ (banana width). The dynamics of the gyro and bounce averaged trapped particles or banana centre distribution function $f_s$ is determined by the kinetic Vlasov equation,

$$ \frac{\partial f_s}{\partial t} - [J_{0,s}\phi, f_s]_{\alpha,\psi} + \frac{\Omega_D(\psi, \kappa)E}{Z_s} \frac{\partial f_s}{\partial \alpha} = 0. \quad (1) $$

where $\phi$ is the electrostatic potential. $Z_s$ is the charge number of the species. $[\cdots]_{\alpha,\psi}$ is the Poisson bracket in the phase space of toroidal precession angle $\alpha$ and poloidal magnetic flux $\psi$ ($\psi \sim -r^2$, is used as a radial co-ordinate, where $r$ is radius of tokamak). The magnetic flux-function is calculated from the integration $\psi(r) \sim -B_{min} \int_0^r \frac{r}{q_0} dr$, where the negative sign stems from the convention that the direction of poloidal magnetic field $B_\theta$ and the normal direction of internal poloidal surface element is opposite to each other. Then $\psi$ is shifted and normalized such that $\hat{\psi}$ is always positive, and $\hat{\psi} \in [0, 1]$ corresponds to a limited radial extent of a tokamak with $\hat{\psi} = 0$ is associated with the edge, towards the last closed flux surface (LCFS), and $\hat{\psi} = 1$ is associated with the core, towards the hottest central region of a tokamak. This model is not valid inside the open field line regions, where most of the trajectories intercept plasma-facing components. The normalization of $\psi$ and other essential physical quantities are given in Tab. I, where dimensionless normalized quantities are noted with hat. However in the main text, the hat notation is omitted for clarity. $E\Omega_D/Z_s = \omega_{d,s}$ is the energy dependent
precession frequency of species s,

$$ \omega_{d,s} = \frac{q(r)}{r} \frac{E}{q_s B_{min} R_0} \bar{\omega}_d, \quad (2) $$

where $q$, and $R_0$ are the safety factor, and major radius of a tokamak, respectively. $q_s$ is the electric charge of the species $s$, $B_{min}$ is the minimal strength of the magnetic field on a field line, $E \equiv \frac{1}{2}m_s v_G^2 + \mu B_G$ is the particle kinetic energy, $\mu \equiv \frac{m_s v_G^2}{2B_G}$ is the magnetic moment, where subscript $G$ refers to the quantities computed at the position of guiding center, and

$$ \bar{\omega}_d = \frac{2\mathcal{E}(\kappa^2)}{\mathcal{K}(\kappa^2)} - 1 + 4s_0(r) \left( \frac{\mathcal{E}(\kappa^2)}{\mathcal{K}(\kappa^2)} + \kappa^2 - 1 \right), \quad (3) $$

where $\kappa = \sqrt{\frac{1-\lambda}{2\epsilon\lambda}}$ is the trapping parameter which can vary between 0 (for deeply trapped particles) to 1 (at separatrix for barely trapped particle). $\epsilon = a/R_0$ is the inverse of the aspect ratio of the tokamak, $a$ is the minor radius of tokamak, and for a tokamak with large aspect ratio usually the quantity $\epsilon < 1$. $\lambda = \mu B_{min}(\psi)/E$ is the pitch angle, $s_0 = \frac{r}{q(r)} \frac{dq}{dr}$ is the magnetic shear, and $\mathcal{K}(\kappa^2)$ and $\mathcal{E}(\kappa^2)$ are the complete elliptic functions of the first and second kind, respectively. The operator $J_{0,s}$ performs two successive averages: the gyro-average and the bounce-average, which are the average over the cyclotron motion and the banana motion, respectively for the species $s$ is (according to “Padé” expression [25]).

$$ J_{0,s} = \left( 1 - \frac{E}{T_{eq,s}(0)} \frac{\delta_{b0,s}^2}{4} \partial_{\psi}^2 \right)^{-1} \left( 1 - \frac{E}{T_{eq,s}(0)} \frac{q^2 \rho_{c0,s}^2}{4L_{\psi}^2} \partial_{\alpha}^2 \right)^{-1}, \quad (4) $$

where $\rho_{c0,s} = \frac{m_s v_G}{q_s B}$ and $\delta_{b0,s} = q\rho_{c0,s}/\sqrt{\epsilon}$ are the Larmor radius and the banana width (in unit of $\psi$) computed at temperature $T_0$. $T_{eq,s}(0)$ is the equilibrium temperature of species $s$ at $\psi = 0$.

Self-consistency is ensured by a quasi neutrality constraint, including a polarization term $\bar{\Delta}_s\phi$, where $\bar{\Delta}$ is a non-isotropic Laplacian operator,

$$ \bar{\Delta}_s = \left( \frac{q\rho_{c0,s}}{L_{\psi}} \right)^2 \frac{\partial^2}{\partial\alpha^2} + \delta_{b,s}^2 \frac{\partial^2}{\partial\psi^2}. \quad (5) $$

where $L_{\psi} = aR_0B_{\theta}$ is the radial length of the simulation box in unit of $\psi$. The quasi-neutrality equation with the approximation $H_{eq,s}(\psi) = E(1 + \Omega_D\psi) \approx E$ reads as follows [25, 30],

$$ \frac{2}{\sqrt{\pi n_{eq}(0)}} \sum_s Z_s \int_0^1 \kappa \mathcal{K}(\kappa^2) d\kappa \int_0^\infty J_{0,s} f_s \sqrt{E} dE \\ = \sum_s \frac{eZ_s^2}{T_{eq,s}(0)} \left[ \frac{1 - f_t}{f_t} (\phi - \epsilon_{\phi,s} \langle \phi \rangle_{\alpha}) - \bar{\Delta}_s \phi \right], \quad (6) $$

<table><caption>TABLE I: Normalization of the plasma parameters. Physical quantities are noted without a hat, and dimensionless quantities are noted with a hat. Here $\omega_{d,0} = q_0T_0/(er_0R_0B_0)$ is a typical precession frequency of strongly trapped ion at $E = T_0$. $n_0$ and $T_0$ are arbitrary normalizing ion density and temperature such that $\hat{n}_s = \hat{T} = 1$ at $\hat{\psi} = 0$. The quantity $L_\psi$ is the radial size of the simulation box in magnetic flux unit. The minor radius $a$, the Larmor radius $\rho_{c0}$, and the banana width $\delta_b$ are all expressed in units of $\psi$. However in the main text, the hat notation is omitted for clarity.</caption><thead><tr><th>Quantity</th><th>e.g.</th><th>Normalization</th></tr></thead><tbody><tr><td>Time</td><td>t, ω<sup>-1</sup></td><td>ħ = ω<sub>d,0</sub>t</td></tr><tr><td>Poloidal magnetic flux</td><td>ψ, a, ρ<sub>c0</sub>, δ<sub>b</sub></td><td>ψ̂(r) = (ψ(r) - ψ(a)) / L<sub>ψ</sub></td></tr><tr><td>Electric potential</td><td>φ</td><td>φ̂ = φ/(ω<sub>d,0</sub>L<sub>ψ</sub>)</td></tr><tr><td>Energy</td><td>E</td><td>Ê = E/T<sub>0</sub></td></tr><tr><td>Density</td><td>n<sub>s</sub></td><td>n̂<sub>s</sub> = n<sub>s</sub>/n<sub>0</sub></td></tr><tr><td>Temperature</td><td>T</td><td>T̂ = T/T<sub>0</sub></td></tr><tr><td>Distribution function</td><td>f<sub>s</sub>, F<sub>eq</sub></td><td>f̂<sub>s</sub> = (1/n<sub>0</sub>) (2πT<sub>0</sub>/m)<sup>3/2</sup> f<sub>s</sub></td></tr></tbody></table>

where, $n_{eq}(0)$ is the equilibrium density at $\psi = 0$, $f_t = \frac{2\sqrt{2\epsilon}}{\pi}$ is the fraction of trapped particles, which scales as $\sqrt{r}$, whereas we consider a constant $f_t$. $\langle \dots \rangle_\alpha$ corresponds to the average on the angle $\alpha$. Here passing particles are treated quasi-adiabatically. In all the previous studies of trapped particle modes instability [8, 25, 26], the gyro-bounce averaged kinetic model was simplified by considering the equilibrium Hamiltonian $H_{eq}(\psi) = E(1 + \Omega_D\psi) \approx E$, which was used to derive the quasi-neutrality equation Eq. (6). This approximation is valid only at $\psi = 0$ and/or $\Omega_D \ll 1$. The exact expression for equilibrium Hamiltonian is,

$$ H_{eq}(\psi, \kappa) = E \left( 1 + \int_{0}^{\psi} \Omega_{D}(\tilde{\psi}, \kappa) d\tilde{\psi} \right) = E \Lambda_{D}(\psi, \kappa), \quad (7) $$

where $\Lambda_D(\psi, \kappa) = 1 + \int \Omega_D(\tilde{\psi}, \kappa) d\tilde{\psi}$, with $\tilde{\psi}$ is a variable of integration. In the previous derivation of gyro-bounce averaged Vlasov equation Eq. (1) the $\psi$ derivative of the equilibrium Hamiltonian was considered as $\frac{dH_{eq}}{d\tilde{\psi}} = E\Omega_D$, which remains unchanged for this new Hamiltonian Eq. (7). Therefore the expression for gyro-bounce averaged Vlasov equation Eq. (1) remains unchanged. With this new equilibrium Hamiltonian Eq. (7) the elementary volume in the phase-space can be written as $d^3v = 4\pi\sqrt{2}m^{-3/2}\sqrt{E}\Lambda_D^{3/2}dE \frac{d\lambda}{4\Omega_D}$. For simplicity we consider, the term $\bar{\Delta}\phi$ associated with the polarization in the quasi-neutrality equation Eq. (6) remain unchanged under this new modification in the equilibrium Hamiltonian. By keeping unchanged the right-hand side of the previous quasi-neutrality equation Eq. (6) and using the modified expression of elementary volume $d^3v$ for the integration of gyro-bounce averaged distribution $f_s$, the
quasi-neutrality equation can be written as

$$ \frac{2\Lambda_D(\psi, \kappa)^{3/2}}{\sqrt{\pi} n_{eq}(0)} \sum_s Z_s \int_0^1 \kappa K(\kappa^2) d\kappa \int_0^\infty J_{0,s} f_s \sqrt{E} dE \\ = \sum_s \frac{eZ_s^2}{T_{eq,s}(0)} \left[ \frac{1-f_t}{f_t}(\phi - \epsilon_{\phi,s}\langle\phi\rangle_\alpha) - \bar{\Delta}_s \phi \right], \quad (8) $$

Therefore there is an additional multiplication term $\Lambda_D(\psi, \kappa)^{3/2}$ that arises in the left-hand side of the previous quasineutrality condition Eq. (6). In this model, the gyro-bounce averaged distribution function $f_s$ bears 4 dimensions $(\alpha, \psi, E, \kappa)$, and in a further reduced limit of a single $\kappa$ value it reduces to a 3-dimension model.

## 1. Global linearised model

Considering an initial very small amplitude perturbation of gyro-bounce averaged distribution $f_s = F_{eq,s} + \tilde{f}_s$ and potential $\phi = \tilde{\phi}$ in Fourier space as $\tilde{f}_s = \sum_{n,\omega} f_{s,n,\omega}(\psi, E, \kappa) \exp\{i(n\alpha - \omega t)\}$ and $\tilde{\phi} = \sum_{n,\omega} \phi_{n,\omega}(\psi) \exp\{i(n\alpha - \omega t)\}$, then neglecting the higher order nonlinear terms $[J_{0,s}\tilde{\phi}, \tilde{f}_s]_{\alpha,\psi}$, the linearised form of Vlasov equation can be written as,

$$ \frac{\partial \tilde{f}_s}{\partial t} - [J_{0,s} \tilde{\phi}, F_{eq,s}]_{\alpha, \psi} + \frac{\Omega_D E}{Z_s} \frac{\partial \tilde{f}_s}{\partial \alpha} = 0. \quad (9) $$

$n$ and $\omega$ are the mode number (along $\alpha$), and angular frequency of the Fourier modes respectively. Here we consider the normalized gyro-bounce averaged equilibrium distribution $F_{eq,s}$ is Maxwellian energy distribution, which is independent of $(\alpha, t)$. Using the new equilibrium Hamiltonian $F_{eq,s}$ can be written as:

$$ F_{eq,s}(\psi, E, \kappa) = \frac{n_{eq,s}(\psi)}{T_{eq,s}^{3/2}(\psi)} \exp\left(-\frac{E\Lambda_D(\psi, \kappa)}{T_{eq,s}(\psi)}\right), \quad (10) $$

where $T_{eq,s}(\psi)$ and $n_{eq,s}(\psi)$ are the temperature and density profiles of the equilibrium distribution function of species $s$, respectively. In this case, the integration of $F_{eq}$ over the velocity space $d^3v$ leads to $n_{eq}(\psi)$. Hereafter we denote $T_{eq}(\psi)$ as $T(\psi)$ and $n_{eq,s}(\psi)$ as $n_s(\psi)$. The term $E\Lambda_D$ allow us to incorporate the radial variation of particle drift velocity. After substituting $\tilde{f}_s$, $\tilde{\phi}$ and $F_{eq,s}$ in Eq. (9), the solution of Vlasov equation in Fourier space become,

$$ f_{n,\omega} = \frac{\kappa_n(\psi) + \kappa_T(\psi) \left( \frac{E\Lambda_D(\psi)}{T(\psi)} - \frac{3}{2} \right) - \frac{E\Lambda_D}{T(\psi)} \kappa_\Lambda(\psi)}{Z_s^{-1} \Omega_D(\psi, \kappa) E - \frac{\omega}{n}} \quad (11) \\ \{J_{0,n,s} \phi_{n,\omega}(\psi)\} F_{eq,s}(\psi, E, \kappa). $$

where $\kappa_n(\psi) = \frac{1}{n_s(\psi)} \frac{dn_s}{d\psi}$, $\kappa_T(\psi) = \frac{1}{T_s(\psi)} \frac{dT_s}{d\psi}$ and $\kappa_\Lambda(\psi) = \frac{1}{\Lambda_D(\psi)} \frac{d\Lambda_D}{d\psi} = \frac{\Omega_D(\psi)}{\Lambda_D(\psi)}$ are the logarithmic gradients of density, temperature and equilibrium Hamiltonian $H_{eq}(\psi)$, respectively.

Considering both electron and ion contributions, the quasi-neutrality condition Eq. (8) can be written as,

$$ \frac{\sqrt{\pi}}{2T_i(0)} [C_{ad}(\phi_{n,\omega} - \epsilon_{\phi}\langle\phi\rangle_{\alpha}) - C_{pol}\bar{\Delta}\phi_{n,\omega}] = \mathcal{N}_{n,i} - \mathcal{N}_{n,e}, \\ \mathcal{N}_{n,i} = \frac{\Lambda_D^{3/2}}{n_i(0)} \int_0^1 \int_0^\infty J_{0,i} f_{i,n,\omega}(\psi, E, \kappa) \sqrt{E} dE d\kappa K(\kappa^2) d\kappa, \quad (12) \\ \mathcal{N}_{n,e} = \tau \frac{\Lambda_D^{3/2}}{n_e(0)} \int_0^1 \int_0^\infty J_{0,e} f_{e,n,\omega}(\psi, E, \kappa) \sqrt{E} dE d\kappa K(\kappa^2) d\kappa, $$

where $\tau = \frac{T_i}{T_e}|_{\psi=0}$, $C_{pol} = \frac{q_i \omega_0 L_\psi}{T_0}$, $C_{ad} = \frac{1-f_t}{f_t}(1-\tau)C_{pol}$, $\epsilon_\phi = \frac{\epsilon_{\phi,i} + \tau\epsilon_{\phi,e}}{1+\tau}$ and $\bar{\Delta}\phi = \bar{\Delta}_i\phi + \tau\bar{\Delta}_e\phi$. In the limit of a constant pitch-angle, the value of $\Omega_D$ is constant along $\kappa$, then the $\kappa$ integration in Eq. (12) can be simplified as, $\int_0^1 \kappa K(\kappa^2) d\kappa = 1$. Substituting the value of $f_{n,\omega}$ from Eq. (11) in Eq. (12), the expression of $\mathcal{N}_s$ becomes

$$ \mathcal{N}_n = \frac{\Lambda_D(\psi)^{3/2}}{n_s(0)} \int_0^\infty \sqrt{E} J_{0,n} \left[ \frac{n_s(\psi)}{T^{3/2}(\psi)} \exp\left(-\frac{E\Lambda_D}{T(\psi)}\right) \frac{\kappa_n(\psi) + \kappa_T(\psi) \left(\frac{E\Lambda_D}{T(\psi)} - \frac{3}{2}\right) - \kappa_\Lambda(\psi) \frac{E\Lambda_D}{T(\psi)}}{Z^{-1}\Omega_D(\psi)(E-\chi)} (J_{0,n}\phi_{n,\omega}) \right] dE, \quad (13) $$

where $\chi = \frac{\omega}{nZ^{-1}\Omega_D(\psi)}$, and $\omega$ has both real and imaginary parts $\omega = \omega_r + i\gamma$. Due to the term $(E - \chi_s)$ in denominator there is a possibility of resonance between wave and particle motion. For ions $Z_i$ is positive, therefore the resonance occurs only when the phase velocity of wave has the same sign as the ion precession drift (ie., $\omega > 0$). By substituting $\bar{\Delta}_s$ from Eq. (5) in the left-hand side of quasineutrality condition Eq. (12) we define the differential operator

$$ C_n = \frac{\sqrt{\pi}}{2T_i(0)} \left[ C_{ad}(1 + \epsilon_{\phi} \delta_{n,0}) + C_{pol} \left\{ (\rho_i^{*2} + \tau \rho_e^{*2})(-n^2) - (\delta_{bi}^2 + \tau \delta_{be}^2) \frac{\partial^2}{\partial \psi^2} \right\} \right], \quad (14) $$

where $\rho_s^* = \frac{q\rho_{c0,s}}{L_\psi}$, and $\delta_{n,0}$ is Kronecker delta with the value 1 for $n=0$ and for $n \neq 0$, it is 0. Therefore the dispersion relation becomes,

$$ C_n \phi_{n,\omega} = \mathcal{N}_{n,i}^* \phi_{n,\omega} - \mathcal{N}_{n,e}^* \phi_{n,\omega} \quad (15) $$

where $\mathcal{N}_{n,s}^* = \frac{\mathcal{N}_{n,s}}{\phi_{n,\omega}}$ is actually a differential operator (Eq. (4) and (13)) acting on $\phi_{n,\omega}$. We will come back to this issue in Sec. III. Hereafter we remove the subscript $\omega$ from $\phi_{n,\omega}$, and denote it as $\phi_n$, because the linear dispersion relates $\omega$ to $n$ values. By substituting the expression of $C_n$ in Eq. (15) one can derive a 2nd order linear differential equation of $\phi_{n,\omega}$ as,

$$ \frac{d^2\phi_n}{d\psi^2} + Q_n(\psi)\phi_n = 0, \\ Q_n(\psi) = \frac{\mathcal{N}_{n,i}^* - \mathcal{N}_{n,e}^* - \frac{\sqrt{\pi}}{2T_i(0)} [C_{ad} + C_{pol}(\rho_i^{*2} + \tau\rho_e^{*2})n^2]^{(16)}}{\frac{\sqrt{\pi}}{2T_i(0)} C_{pol}(\delta_{bi}^2 + \tau\delta_{be}^2)}. $$
Considering $T_i(0) = T_0$, $n_i(0) = n_e(0) = n_0$, in normalized unit those become $T_i(0) = n_s(0) = 1$. Here we study the modes for which $n \neq 0$, therefore the term $\epsilon_\phi \delta_{n,0}$ in Eq. (14) vanishes. For simplicity to study specifically the trapped ion mode instability (TIM) we will neglect the electron perturbation $\hat{f}_e = 0$ which leads to $\mathcal{N}_e^* = 0$. However, the same method is applicable to trapped electron mode (TEM) instability by setting $\mathcal{N}_i^* = 0$ and $\mathcal{N}_e^* \neq 0$. Moreover, for $\tau \le 1$, $\rho_i^{*2} \gg \rho_e^{*2}$ and $\delta_{bi}^2 \gg \delta_{be}^2$. Therefore after neglecting those terms associated with electrons contribution and substituting $T_{eq,i} = 1$ in eq. (16), the modified expression for $Q_n(\psi)$ becomes,

$$Q_n(\psi) = \frac{\mathcal{N}_{n,i}^* - \frac{\sqrt{\pi}}{2} \left[ C_{ad} + C_{pol} \rho_i^{*2} n^2 \right]}{\frac{\sqrt{\pi}}{2} C_{pol} \delta_{bi}^2} \quad (17)$$

The contributions of the density gradient, temperature gradient and the gradient in the equilibrium Hamiltonian is contained in the expression of $\mathcal{N}_n^*$. Due to the term $-\kappa_\Lambda$, while the increase in logarithmic gradients of temperature and density $\kappa_T$ and $\kappa_n$ help to enhance the TIM instability, the logarithmic gradient of equilibrium Hamiltonian $\kappa_\Lambda$ helps to stabilize the TIM instability. In the limit $H_{eq} \approx E$, the term $\kappa_\Lambda = 0$ and $\Lambda_D = 1$, in Eq. (13). Therefore the growth-rate $\gamma$ of all the modes $n$ of TIM instability for the new modified equilibrium Hamiltonian $H_{eq}(\psi, \kappa, E)$, will be significantly smaller compared to the case with the limit $H_{eq} \approx E$. The dispersion relation of the TIM mode instability in the limit of $H_{eq} \approx E$ is presented in the Appendix A. Since the solution in Fourier space Eq. (11) $f_{n,\omega} = 0$ for $n = 0$, the linear analysis is valid only for the mode numbers $n > 0$. However, $n = 0$ mode is linearly stable, it cannot extract free energy from the equilibrium gradients.

## 2. Local linearised analysis

The local linear stability analysis of this reduced gyro-bounce averaged model with the new modified equilibrium Hamiltonian can be obtained by expanding $F_{eq,s}(\psi, \kappa, E)$ eq. (10) up to 1st order of $\psi$ around $\psi = 0$, and substituting $\partial_\psi^2 \phi = -k^2 \phi$ (with $k = \pi$) in the dispersion relation, which leads to the new simplified dispersion relation as, $C_n - \mathcal{N}_{n,i}^* = 0$. From this, using Plemelj formula [31], the threshold frequency value of the real part of frequency $\omega_r$ for TIM instability can be derived as:

$$\omega_r^{\text{new}} = \left( \frac{\frac{3}{2} \kappa_{T0} - \kappa_{n0}}{\kappa_{T0} - \kappa_{\Lambda 0}} \right) \frac{\Omega_{D0}}{\Lambda_D} T_{0n}, \quad (18)$$

where the subscript '0' denotes the variables value at $\psi = 0$. Here we restrict the analysis to the case of resonant interactions only, i.e. $\omega_r > 0$, hence Eq. (18) is valid for the cases with $(3\kappa_T/2 - \kappa_n)/(\kappa_T - \kappa_\Lambda) > 0$. The threshold value of $\kappa_T$ for the instability can be written

as

$$\kappa_{T,th}^{\text{new}} = \frac{C_n \Omega_{D0}}{\Lambda_{D0} \int_{0}^{\infty} J_{0,n}^2 \exp(-\xi) \sqrt{\xi} d\xi} + \kappa_{\Lambda 0} \quad (19)$$

where $\xi = \frac{E\Lambda_D}{T}$. The threshold values of $\omega_r$ and $\kappa_T$ Eq. (A5-A6) in the limit of $H_{eq} \approx E$ [9] can be recovered by substituting $\Lambda_D = 1$ and $\kappa_\Lambda = 0$ in Eq. (18-19). The relation between these two threshold $\kappa_T$ values is

$$\kappa_{T,th}^{\text{new}} = \frac{\kappa_{T,th}}{\Lambda_{D0}} + \kappa_{\Lambda 0}, \quad (20)$$

where $\kappa_{T,th}^{new}$ and $\kappa_{T,th}$ are the threshold values of $\kappa_T$ for TIM instability in the case with $\kappa_\Lambda$ and without $\kappa_\Lambda$, respectively. For $\Lambda_{D0} \approx 1$, $\kappa_{T,th}^{new} \gg \kappa_{T,th}$, which suggests a relatively strong gradient in temperature profile is required in order to obtain TIM instability, compared to the case with $H_{eq} \approx E$. In case of adiabatic electron response there is a non-resonant branch of TIM for $\chi < 0$ [25]. But here we are focusing only on the resonant branch of TIM, for which $\chi > 0$.

In the present study, unlike the local linear analysis, instead of Taylor expansion of $F_{eq,s}$ around $\psi = 0$, a full $\psi$ dependent Maxwellian distribution for the equilibrium gyro-bounce averaged particle distribution $F_{eq,a}(\psi, \kappa, E)$ is considered. Therefore using this global linear analysis, effects of any type of temperature, density and precession frequency profiles on the stability of trapped ion modes (TIM) can be investigated. In the next section we will solve the differential equation of $\phi$ eq. (16) with the expression for $Q_n$ from Eq. (17) using spectral method.

## III. GLOBAL LINEAR ANALYSIS: SPECTRAL METHOD

The $Q_n$ function in Eq. (17) depends on $\psi$ in a intricate manner. However in the limit of local linear analysis [9], $Q_n$ is independent of $\psi$. Therefore $-Q_n$ becomes the eigenvalue of the operator $\partial_\psi^2$ in Eq. (16). But in our global linear analysis, due to full Maxwellian distribution $F_{eq}$, it is impossible to eliminate the $\psi$ dependency in $Q_n$, ie., $-Q_n$ is no more an eigenvalue of $\partial_\psi^2$. Here we use a spectral approach [28, 29, 32] to solve Eq. (16). We can consider our co-ordinate system $(\psi, \alpha)$ as a polar co-ordinate system with $\psi$ as radial axis and $\alpha$ as angular axis. If we consider the solution $\phi(\psi, \alpha)$ as a scalar analytic function in the region $0 \le \psi \le 1$, according to the Theorem-1 of the page-374 in [32], using a spectral approach the solution $\phi(\psi, \alpha)$ can be expanded as:

$$\phi(\psi, \alpha) = \sum_{n=-\infty}^{\infty} \psi^{|n|} \sum_{m=1}^{\infty} \hat{\phi}_{n,m} \mathfrak{F}_m(\psi) \exp(in\alpha) \quad (21)$$

where $\mathfrak{F}_m(\psi)$ is an even power function of $\psi$, and the sum index $m$ has values $m = 1, 2, \dots, \infty$. In practice, a truncation is performed so that it is limited to the first $M$ number functions of $\mathfrak{F}_m$. In Eq. (16) all the differentiation in $\alpha$ is taken in Fourier-space. Therefore according
to the spectral approach from Eq. (21) the solution of Eq. (16) $\phi_n(\psi)$ can be expressed as:

$$ \phi_n(\psi) = \psi^{|n|} \sum_{m=1}^{M} \hat{\phi}_{n,m} \mathfrak{F}_m(\psi) \quad (22) $$

The boundary conditions for $\phi_n(\psi)$ in Eq. (16) are $\phi|_{\psi=0} = 0$ and $\phi|_{\psi=1} = 0$. Our first choice for $\mathfrak{F}_m$ was $\mathfrak{F}_m = (1 - \psi^2)\psi^{2(m-1)}$, but this choice gives poorly conditioned matrices. For good conditions matrices the coefficients $\hat{\phi}_{n,m}$ go down towards zero exponentially with increasing $m$ and at $m = M$, $\hat{\phi}_{n,M} \approx 0$. Then we consider orthogonal polynomials, Chebyshev polynomials of first kind $T_{2(m-1)}(\psi)$, for constructing $\mathfrak{F}_m$ as:

$$ \mathfrak{F}_m(\psi) = (1 - \psi^2)T_{2(m-1)}(\psi). \quad (23) $$

which were used for studying the collisional drift wave and ITG instabilities in a cylindrical plasma, using spectral method [28, 29]. However, in case of TIM instability we found that this choice of function gives good conditions matrices for lower values of $n$. For higher values of $n$, the coefficients $\hat{\phi}_{n,m}$ do not go down towards zero exponentially with increasing $m$. Then we choose $n = 1$ for the factor $\psi^{|n|}$, in the expansion function Eq. (22) and the dependency of solution on the mode-number $n$ is entirely determined by the expression $Q_n$ in Eq. (16). With this new choice the coefficients $\hat{\phi}_{n,m}$ go down towards zero exponentially with increasing $m$. Moreover we have checked that, using this new choice of expansion function, one can recover the solutions of collisional drift wave instability study, as given in the references [28, 29].

By defining $C_m = \psi \mathfrak{F}_m$ the expression of $\phi_n(\psi)$ can be written as,

$$ \phi_n(\psi) = \sum_{m=1}^{M} \hat{\phi}_{n,m} C_m. \quad (24) $$

Here $C_m = \psi(1-\psi^2)T_{2(m-1)}(\psi)$. Using this new choice of $\phi_n(\psi)$ rather good conditioning properties of matrices are obtained for large values of $n$ in our study. Substituting the value of $\phi$ in Eq. (16) we get,

$$ \sum_{m=1}^{M} \hat{\phi}_{n,m} \frac{d^2 C_m(\psi)}{d\psi^2} = - \sum_{m=1}^{M} \hat{\phi}_{n,m} Q_n(\psi) C_m(\psi), \quad (25) $$

where $m$ is the index of the spectral function $C_m$, and $n$ is wave number of TIM instability along the toroidal precession angle $\alpha$. The 2$^{nd}$-order differential Eq. (25), is solved numerically. By proceeding this way, Eq. (25) is only solved at the points $\psi$, where its left hand side is vanishing for the first dropped higher order term of the expansion, namely the term $C_{M+1}(\psi) = \psi(1-\psi^2)T_{2M}(\psi)$. Those $\psi$ values are called the collocation points, which are defined by the $M$ zeros of the Chebyshev polynomial $T_{2M}(\psi)$ situated in the interval $0 \le \psi \le 1$, ie.,

$\psi_l = \cos\left(\frac{(2l-1)\pi}{4M}\right)$ for $l \in \{1, \dots, M\}$. Since $C_{M+1}$ vanishes at these $\psi$ locations, for $J_0 = 1$ the right hand side of Eq. (25) also vanishes at $m = M + 1$. Moreover, for gyro-bounce average operator $J_0$ according to Eq. (4) we solve the Eq. (25) iteratively where the operator $N_n^*$ inside $Q_n(\psi)$ is applied on a solution $\phi_{n-\Delta n}(\psi)$ which is obtained from the previous iteration. Therefore, at any iteration $Q_n(\psi)$ is a known function of $\psi$, it does not directly operate on $C_m$. Hence at the collocation points the right hand side of Eq. (25) vanishes for the $m = (M+1)$. For minimizing the truncation error, the discretization in $\psi$ is considered at those collocation points. We will explain this iteration method at the end of this section. For each value of wave number $n$ along $\alpha$, Eq. (25) is evaluated at these collocation points and yields the matrix problem,

$$ \mathcal{M}_D \hat{\phi}_n = -\mathcal{M}_Q \hat{\phi}_n \quad (26) $$

where $\hat{\phi}_n = \begin{pmatrix} \hat{\phi}_{n,1} \\ \vdots \\ \hat{\phi}_{n,M} \end{pmatrix}$ is the vector representing the solution, and the matrices are,

$$ (\mathcal{M}_D)_{l,m} = \frac{d^2}{d\psi_l^2} C_m(\psi_l), \quad (27) $$

$$ (\mathcal{M}_Q)_{l,m} = Q_n(\psi_l)C_m(\psi_l), \quad (28) $$

with $l, m \in \{1, \dots, M\}$. Then we scan the $(\omega_r, \gamma)$ plane (where $\omega = \omega_r + i\gamma$), and search for values of $\omega$ such that one eigenvalue of the matrix $M = \mathcal{M}_D + \mathcal{M}_Q$ vanishes within machine precision. During this search, there is no restriction on the range of $\omega_r$ values, and we search only for positive $\gamma$ values. For this purpose a method that finds the minimum of a scalar function of several variables, starting at an initial guess value, and iterates using the *simplex search method* [33] is used. Then using the eigenvector of matrix $M$ associated with this smallest eigenvalue, the electric potential $\phi_n(\psi)$ can be constructed using Eq. (24). The solution $\phi_n(\psi)$ has both real and imaginary part. The *spectral-convergence* or rate of change of the coefficients $\hat{\phi}_{n,m}$ values with different $M$ values depends on the nature of the matrix $\mathcal{M}_Q$. Usually acceptable smooth $\hat{\phi}_{n,m}$s follows exponentially decaying functions with $m$, i.e., $\hat{\phi}_{n,m} = \hat{\phi}_{n,0} \exp(-|\beta|m)$, where $\beta$ is a constant. For different initial guess values the minimum searching method can converge towards different values of $\omega$. Among these solutions one finds the couple $(\omega_r, \gamma)$ for which the instability growth rate $\omega_i$ is maximum. We will take that particular set as a solution $(\omega_r, \gamma)$. For generating the curve $(\omega_r$ Vs. $n)$ and $(\gamma$ Vs. $n)$ we start from mode $n = 1$ and search for the highest value of growth rate, by choosing different initial guess values as $\omega_{r0} \in [0, 3\omega_{rL}]$ and $\gamma_0 \in [0, \omega_{rL}]$ (since $\omega_r \ge \gamma$), and with the intervals $\Delta\omega_{r0} = 0.1\omega_{rL}$, $\Delta\gamma_0 = 0.01\omega_{rL}$, where $\omega_{rL}$ is the threshold value of $\omega_r$ obtained from the local linear analysis Eq. (18). Actually the range and the
intervals for the initial guess values depend on the nature of instability ie., on the matrix $M_Q$. After finding the highest value of growth rate for $n = 1$, increase the $n$ value slowly with a step $\Delta n = 0.1$. Though the mode numbers $n$ along $\alpha$, are integer number, there is no such restriction on $n$ for solving the differential Eq. (25) for fractional values of $n$. In the search method, we choose the solution $\omega$ of $n - \Delta n$ step, as initial guess $(\omega_{r0}, \gamma_0)$ value to search the solution for mode number $n$. Here we consider the $(\omega_r, \text{Vs. } n)$ and $(\gamma, \text{Vs. } n)$ profiles vary smoothly with $n$, such that small change in $n$, changes the $M$ matrix by small amount, and helps to find-out the values $(\omega_r, \gamma)$ close to the initial guess value $(\omega_{r0}, \gamma_0)$, and keep the solutions within the same branch with highest growth rate. Larger step in $\Delta n$ in the searching method, may lead to a departure of the solution $(\omega_r, \omega_i)$ from the branch with highest growth rate. However if the solution departs from any particular branch, i.e. if it leads to a sudden change in potential profile for two consecutive $n$ values, then one has to decrease the $\Delta n$ value.

The main difficulty for solving the Eq. (16) using the spectral method is to construct the $Q_n(\psi)$ matrix. Because for calculating $N_{n,i}^*$, the gyro-bounce average operator $J_{0,n,i}$ is applied to the potential solution $\phi_n$ in Eq. (13), which is an unknown quantity of Eq. (16). Therefore we will solve the Eq. (16) for $\phi$ iteratively. Moreover to calculate the gyro-bounce average of any function $\mathcal{F}(\psi)$ using the Padé expression Eq. (4), we solve the differential equation,

$$ \left( 1 - \frac{E}{T(0)} \frac{\delta_{b0,s}^2}{4} \frac{d^2}{d\psi^2} \right) \bar{\mathcal{F}}(\psi) = \left( 1 + \frac{E}{T(0)} \frac{\rho^{*2}}{4} n^2 \right)^{-1} \mathcal{F}(\psi) \quad (29) $$

where, $\mathcal{F}$ is a known function, on which the gyro-bounce average operator $J_0$ is applied and generates the gyro-bounce averaged quantity $\bar{\mathcal{F}}$. We solve the Eq. (29), using the spectral method with the boundary condition $\bar{\mathcal{F}} = 0$ at $\psi = 0$ and $\psi = L_\psi = 1$. Finally Eq. (16) is solved iteratively, by calculating $N_n$ for $\phi_n$ solution from the previous iteration. For the 1st iteration with mode number $n$, the potential solution $\phi_{(n-\Delta n)}(\psi)$ of the mode $n - \Delta n$ is taken into account to calculate $N_n$. But for the calculation with mode $n = 1$, at the beginning (1st iteration) the potential solution is not available, therefore at first, Eq. (16) is solved using the spectral method for the case with $J_{0,s} = 1$ and construct the potential using Eq. (24). Then the gyro-bounce average operator Eq. (4) is applied on that potential solution for calculating $N_n$ at the 1st iteration for $n = 1$, and construct the solution $\phi_n$. In the next iteration this $\phi_n$ is used to calculate $N_n$ and update the solution $\phi_n$. After every iteration, $\Delta\phi_n$, the difference in solution $\phi_n$ from its value at previous iteration, is calculated and this iteration method continuous until achieving the accuracy limit $|\Delta\phi_n| \sim 10^{-3}$.

Though this solver based on spectral method is quite robust for studying the linear stability of TIM for different dependent parameters (temperature, density and precession frequency) profiles, it has some limitations.

We have observed this spectral method can calculate the $(\gamma, \omega_r)$ precisely for temperature profiles having $\kappa_T < 10\kappa_{T,th}$, and beyond that the method becomes unstable, and even a very small change in $\Delta n < 0.001$ departs the solution from the expected branch with highest growth-rate. $\kappa_{T,th}$ is the threshold value of instability which is given by Eq. (19 and A6). For $H_{eq} \approx E$ and $\kappa_n = 0$ the threshold value is $\kappa_{T,th} \sim 0.12$ for the mode number $n = 1$ [9]. Moreover, this method can be applied only for the potential solutions having zero values at the boundary $(\phi = 0$ at $\psi = 0$ and $\psi = 1)$. In the next two sections Sec. IV-V, this spectral technique is used for studying the effect of temperature and precession frequency profiles (having $\kappa_{T,max} < 10\kappa_{T,th}$) on TIM instability for the case in the limit $H_{eq} \approx E$ and for the exact equilibrium Hamiltonian $H_{eq}(\psi, \kappa, E)$. $\kappa_{T,max}$ is the maximum value of a $\kappa_T(\psi)$ profile.

## IV. EFFECT OF TEMPERATURE AND PRECESSION FREQUENCY ON TIM INSTABILITY IN THE LIMIT $H_{eq} \approx E$

In this section we consider the case with the limit $H_{eq} \approx E$ and present the effect of different temperature profiles, and different precessional frequency profiles $\Omega_D(\psi, \kappa)$, corresponding to different safety factor profiles, on the trapped ion mode instability. We can substitute $\Lambda_D = 1$ and $\kappa_\Lambda = 0$ in the Eq. (13) and get the required expression for $N_n$ Eq. (A3), in the limit $H_{eq} \approx E$. In this limit the growth-rate mostly depends on the terms $\kappa_T$ and $\kappa_n$ in Eq. (13). Moreover, for all the cases the results from the global linear analysis is compared with the linear TERESA simulations. Throughout this paper we have considered a flat density profile $\kappa_n = 0$, which prevents the generation of electron roots of trapped ion mode (TIM) $(\omega_r < 0)$ with propagation along electron diamagnetic direction (appendix of [34]). Here all the modes are ion roots of TIM which propagate along the ion diamagnetic direction. In this sense the modes are “pure-TIM”.

### A. Different temperature profiles

In this sec. (IV A) we use four different types of temperature profiles and a flat density profile $\kappa_n = 0$, and investigate their effect on the trapped ion mode (TIM)
instability. These temperature profiles are presented as,

$$
\begin{align}
T_1(\psi) &= T_0 (1 + G_T \psi) \nonumber \\
T_2(\psi) &= T_0 + \frac{G_T}{2} \left( 1 + L_1 \left[ \log \left( \cosh \frac{\psi - \psi_1}{L_1} \right) \right. \right. \nonumber \\
&\qquad \left. \left. - \log \left( \cosh \frac{\psi - (1 - \psi_1)}{L_1} \right) \right] \right) \tag{30} \\
T_3(\psi) &= T_0 \exp(G_T \psi) \nonumber \\
T_4(\psi) &= T_0 \exp \left[ G_T R_T \tanh \frac{\psi^2 - \psi_1^2}{R_T^2} \right] \nonumber
\end{align}
$$

Here we use the constant parameters $T_0 = 1$, $G_T = 0.25$, $L_1 = 0.025$, $\psi_1 = 0.2$, $R_T = 0.5$ and for the temperature profile $T_4$ we have chosen two different values of $\psi_1$, (0.2 and 0.5). Fig. 1 presents all these five different temperature profiles. $T_1(\psi)$ is the temperature profile with constant temperature gradient. $T_2(\psi)$ has constant gradient inside the simulation domain from $0.2 \le \psi \le 0.8$. The width of this region is controlled by the parameter $\psi_1$ in $T_2$. Near the boundary, the gradient is zero along the $\psi$ direction. Therefore there is a sudden change in gradient at $\psi = 0.2$ and 0.8. This type of profile with zero gradient at boundary is used in TERESA simulation to prevent numerical instabilities arising from the boundaries. $T_3(\psi)$ is an exponentially increasing temperature profile. $T_4(\psi)$ is a special type of temperature profile which has zero gradient near boundary, and inside the box it has finite gradient which can be controlled by the parameter $\psi_1$ keeping $G_T$ constant. Here two different values of $\psi_1$ are used $\psi_1 = 0.2$ and $\psi_1 = 0.5$. For both values of $\psi_1$ the gradient changes smoothly from zero to finite value as $\psi$ goes from boundary towards center region. Therefore this profile has less-chance to generate numerical instabilities. Temperature $T_i$ is normalized to $T_0$, and the poloidal magnetic flux is normalized such that $\psi = 0$ is close to last-closed flux surface and $\psi = 1$ is in the core region of a tokamak. Therefore $T_i = 0.9$ is associated with the ion temperature $0.9T_0$. Moreover, for any particular temperature profile, as an example $T_1$, the core-most temperature ($\psi = 1$) is $1.25T_0$ and the temperature at edge ($\psi = 0$) is $T_0$, so that $\psi \in [0, 1]$ corresponds to a limited radial extent of a tokamak where the temperature changes by only 25%. In tokamak experiments the ion temperature profile depends on operation-mode, during L-mode operation it looks like a bell shaped profile with a flat top near the core and decreases towards the Scrape-Off Layer region, and during H-mode operation due to generation of an edge transport barrier the plasma temperature profile changes significantly. It is difficult to present the entire temperature profile using a single function, but the profiles $T_1$, $T_2$, $T_3$ and $T_4$ are somewhat representative of radial sections as we have checked from different tokamak experiments e.g. data from JET [3, 35], COMPASS [36].

Fig. 2 presents the $\kappa_T(\psi) = \frac{1}{T(\psi)} \frac{dT}{d\psi}$ profiles for all these five temperature profiles. The $\kappa_T$ profiles change

<figure><img src="image_1.png" /><figcaption>FIG. 1: Five different types of temperature profiles, $T_1(\psi)$ (Red solid line), $T_2(\psi)$ (black dashed line), $T_3(\psi)$ (blue solid line), $T_4(\psi)$ with $\psi_1 = 0.2$ (green solid line) and $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan solid line). The radial interval in between $\psi = 0$ and $\psi = 1$ covers a certain radial domain in between the last closed fluxed surface (LCFS) and the very core region of a tokamak, respectively. $\rightarrow$ near $\psi = 1$: indicates direction towards core region of tokamak. $\leftarrow$ near $\psi = 0$: indicates direction towards edge region (last closed flux surface) of a tokamak.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>FIG. 2: $\kappa_T = \frac{d \log T(\psi)}{d\psi}$ profiles for the five different temperature profiles.</figcaption></figure>

significantly for all these five temperature profiles. The temperature profile $T_4$ with $\psi_1 = 0.5$ and $0.2$ have $\kappa_T$ profiles with peak values $0.53$ at $\psi \approx 0.6$ and $0.31$ at $\psi \approx 0.4$, respectively and then gradually decreases to zero towards the boundaries. The exponential temperature profile $T_3$ has constant $\kappa_T$. Profiles $T_1$ and $T_2$ have exactly similar $\kappa_T$ profiles, $\kappa_T \propto \psi^{-1}$ within the region $0.2 \le \psi \le 0.8$, and near the boundaries $\psi < 0.2$ and $\psi > 0.8$, $\kappa_T$ suddenly jumps to zero for $T_2$. Therefore in $T_2$ the TIM instability arises due to the temperature profile within $0.2 \le \psi \le 0.8$.
<table><caption>TABLE II: Main input parameters used for studies in sec. IV A.</caption><thead><tr><th>δ<sub>b,i</sub></th><th>ρ<sub>c,i</sub></th><th>C<sub>ad</sub></th><th>C<sub>pol</sub></th><th>Ω<sub>D</sub></th><th>T<sub>eq</sub></th><th>n<sub>eq</sub></th><th>M</th><th>G<sub>T</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.03</td><td>0.1</td><td>0.1</td><td>1</td><td>1</td><td>1</td><td>50</td><td>0.25</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>FIG. 3: Growth rate $\gamma$ Vs. mode number $n$ profiles with $J_0 = 1$, for the five different temperature profiles, $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan line), $T_4(\psi)$ with $\psi_1 = 0.2$ (green line), $T_3(\psi)$ (blue line), $T_1(\psi)$ (red line) and $T_2(\psi)$ (magenta line). The dots denotes the results from the linear TERESA simulation with $J_0 = 1$. Each set of dots marked by same color is from a single simulation.</figcaption></figure>

1. $J_0 = 1$ Case

Since the case with gyro-bounce average operator $J_0 = 1$ is numerically the simplest case, we first validate our global linear analysis for the case with $J_0 = 1$ and with five different temperature profiles. Moreover for mode number $n = 1$, the potential solution from the $J_0 = 1$ case will be used to calculate $N_n$. The other main parameters which are used in this study are presented in the table II. The growth rate $\gamma$ of TIM instability for different mode numbers $n$ is presented in Fig. 3 for all the five temperature profiles. This spectral method is also able to find out the negative growth-rates for the higher mode numbers. Since we are interested only in the unstable TIM modes with positive growth-rates ($\gamma > 0$), in all the figures of $\gamma - n$ profiles, only the modes with $\gamma > 0$ are presented. For the temperature profile $T_4$ with $\psi_1 = 0.5$, the effective temperature gradient and hence the parameter $\kappa_T$ has the highest value compared to other temperature profiles. Therefore as expected the growth-rate of the TIM instability is higher compared to the other cases. For the other temperature profiles the growth-rate decreases as the effective value of $\kappa_T$ decreases. The growth-rate of the highest growing mode for the temperature profiles $T_4$ with $\psi_1 = 0.5$, $T_4$ with $\psi_1 = 0.2$, $T_3$, $T_2$ and $T_1$ are 22.5, 10, 8.3, 5.6 and $6.6\omega_0$, respectively, and the location of most unstable modes are $n \approx 27, 20, 17, 15$ and $16$, respectively. Therefore even

<figure><img src="image_4.png" /><figcaption>FIG. 4: Potential solutions $|\phi_n(\psi)|/\phi_{max}$ for mode $n = 1$ with $J_0 = 1$ for the temperature profiles $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan line), $T_4(\psi)$ with $\psi_1 = 0.2$ (green line), $T_1(\psi)$ (red line), $T_2(\psi)$ (magenta line), and $T_3(\psi)$ (blue line).</figcaption></figure>

for the gyro average operator $J_0 = 1$ case the growth rate of the TIM instability strongly depends on the temperature profile, and the mode with highest growth rate shifts towards zero with decreasing value of $\kappa_T$. For low mode numbers the growth-rate increases with increasing $n$. At very high $n$ (after the mode with highest growth rate) the growth-rate decreases and goes down to zero due to the presence of polarization drift term within the model, and finally generates a bell shaped growth-rate profile. Moreover, the presence of gyro-average operator $J_0$ according to Eq. (4), will enhance this effect that decreases the growth-rate further for higher mode numbers, which is presented in the next sec. IV A 2. All the growth-rate profiles at large values of mode-number $n$, where $\gamma \to 0$ and $\omega_r$ has very large values ($\omega_r \gg \gamma$), slightly depart from expected trend. The reason is that, the integration of Eq. (13) for calculating $N_n$ is a contour integration, and in the complex plane the contour/path of the integration depends on the sign of imaginary part of $\chi$ ie., $\gamma$. Since we are interested in the modes with $\gamma > 0$, we have considered the integration path for the pole at $E = C(\omega_r + i\gamma)$, where $C = (n\Omega_D)^{-1}$. For negative values of $\gamma$ and $\gamma = 0$ the contours of the integration will be different, because the poles are situated at $E = C(\omega_r - i\gamma)$ and $E = C\omega_r$, respectively. One has to take into account these new contours for getting the exact values of growthrate when $\gamma \to 0$ and $\gamma \le 0$. However for higher $\omega_r$ the phase velocity of the wave is very high, and the kinetic effect of trapped ions to TIM modes are negligible. Therefore the growthrate of those modes is left out of our present study. The dots on the growth-rate curves present the measured growth-rates of the different modes $n$ for different temperature profiles from the linear TERESA simulation with $J_0 = 1$. These good agreements of the growth-rate profiles validate the accuracy of our spectral method to the semi-Lagrangian based Vlasov simulation results. For generating an entire
(γ Vs. n) curve, our spectral method based solver takes only 10 – 15 minutes for J₀ = 1 and 20 – 30 minutes for J₀ according to Padé expression, whereas a serial version of TERESA simulation takes around 16 hours. Fig. 4 presents the potential solution φₙ(ψ) profiles for n = 1 for the different temperature profiles. Depending on the κT profiles different φₙ(ψ) profiles are generated. Since the potential solution has both real and imaginary parts, for all the potential profiles within this paper, we take its absolute values and then normalize by its maximum value. All the results that are presented hereafter use the gyro-bounce average operator according to the Padé expression Eq. (4)

## 2. J₀ according to Padé expression

After validating the $J_0 = 1$ case, we intend to solve the differential equation Eq. (16) by calculating $N_n$ from Eq. (A3) with the Padé expression for the gyro-bounce average operator $J_0$ Eq. (4). The gyro-bounce averaged quantity ($J_0f_n$ and $J_0\phi_n$) is calculated by solving the differential Eq. (29) using the spectral method. First $J_{0,n}\phi_{n,\omega}$ is calculated using the solution $\phi_n$, calculated from Eq. (16) by implementing $J_0 = 1$. Then we use this gyro-bounce averaged potential $\bar{\phi} = J_{0,n}\phi_{n,\omega}$ in Eq. (A2) to construct the particle distribution $f_n$, and finally the gyro-bounce averaged distribution function $\bar{f}_n$ is obtained similarly by solving the differential Eq. (29) which is used to calculate $N_n$. Fig. 5 presents the growth-rate of different modes $n$ of TIM instability, for the five different temperatures profiles Eq. (30). In comparison with the previous $J_0 = 1$ cases, the gyro-bounce averaged operator decreases the growth rate $\gamma$ of the TIM instability by a significant amount for all the temperature profiles, roughly by a factor two in terms of highest growth rate. Also the modes with the highest growth rate for each temperature profile are different than the previous $J_0 = 1$ case and they are shifted towards lower mode numbers. Therefore the gyro-average operator reduces the instability of the TIM modes compared to $J_0 = 1$ case. But the effect of the variation in temperature profiles on their growth rate remains unchanged. Similar to previous $J_0 = 1$ case, the growth rate of TIM modes depends on $\kappa_T$ value, and the profile $T_4$ with $\psi_1 = 0.5$ generates the highest growth-rate, which decreases as $\kappa_T$ value decreases for the other temperature profiles. Each dots with different colors on the solid lines present the growth-rate of different modes $n$ for different temperature profiles from the linear TERESA simulation using Padé expression for gyro-bounce average operator. Fig. 6 presents the frequency $\omega_r$ (real part) of different trapped ion modes (TIMs) $n$ for five different temperature profiles. For the lower mode numbers $n \le 10$ frequencies $\omega_r$ are almost similar for all the temperature profiles. The black dashed line presents $\omega_r = \frac{3}{2}n$ which is the solution of the local linear analysis Eq. (A5) in the limit of $\psi \to 0$ [9]. Therefore the frequencies $\omega_r$ of TIMs are higher com-

<figure><img src="image_5.png" /><figcaption>FIG. 5: γ Vs. n profiles with $J_0$ Eq. (4) for the five different temperature profiles, $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan line), and $T_4(\psi)$ with $\psi_1 = 0.2$ (green line), $T_3(\psi)$ (blue line), $T_1(\psi)$ (red line) and $T_2(\psi)$ (magenta line). The dots denote the results from the linear TERESA simulation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 6: Real part of TIM frequency $\omega_r$ Vs. mode-number $n$ with $J_0$ Eq. (4), for five different temperature profiles, $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan line), $T_4(\psi)$ with $\psi_1 = 0.2$ (dashed green line), $T_2(\psi)$ (black line), $T_1(\psi)$ (red line) and $T_3(\psi)$ (blue line). The red+ marker presents the linear TERESA simulation results for $T_4$ with $\psi_1 = 0.5$. The dashed black line presents $\omega_r = 3/2n$, and the black dotted line presents $\omega_r = 2n$.</figcaption></figure>

pared to the local linear case. Indeed we found for these cases they follow $\omega_r \sim 2n$ relation for the global linear analysis. The red marker '+' presents the results from the linear TERESA simulation for the temperature profile $T_4$ with $\psi_1 = 0.5$. Fig. 7 presents the potential solutions $\phi_n(\psi)$ of the Eq. (16) for $n = 1$, for all the five different temperature profiles. The potential profiles are slightly different from the previous case with $J_0 = 1$ (Fig. 4), due to the gyro-bounce average operator $J_0$ Eq. (4). Fig. 8 shows the $\phi_n(\psi)$ solution profile Eq. (24) of mode number $n = 15$ for the temperature $T_4$ with $\psi_1 = 0.5$, which is almost similar to the mode number $n = 1$ (cyan
<figure><img src="image_7.png" /><figcaption>FIG. 7: Potential solutions $|\phi_n(\psi)|/\phi_{max}$ for mode $n = 1$ with $J_0$ Eq. (4) for the five different temperature profiles $T_4(\psi)$ with $\psi_1 = 0.5$ (cyan line), $T_4(\psi)$ with $\psi_1 = 0.2$ (green line), $T_1(\psi)$ (red line), $T_2(\psi)$ (magenta line), and $T_3(\psi)$ (blue line).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 8: Potential solution for the temperature profile $T_4$ with $\psi_1 = 0.5$ for the highest growing mode $n = 15$. Black dashed line presents the potential solution from the TERESA linear simulation.</figcaption></figure>

line in fig. 7). Dashed black line presents the potential profile from the linear TERESA simulation for the temperature profile $T_4$ with $\psi_1 = 0.5$ and mode number $n = 15$. The potential solution is also in good agreement with the TERESA simulation. Fig. 9 presents the decimal logarithms of the coefficients $\hat{\phi}_{n,m}$ of the function $C_m$ in Eq. (24) for the mode number $n = 15$ for the case with temperature profile $T_4(\psi_1 = 0.5)$. The coefficients $\hat{\phi}_{n,m}$ decrease as $m$ increases. This confirms the spectral convergence of $\hat{\phi}_{n,m}$, which follows $|\hat{\phi}_{n,m}| \sim A \exp(-\beta m)$, with $\beta = 0.7$ for $1 \le m \le 18$ (red dashed line in fig. 9). Note that $|\hat{\phi}_{n,m}/\hat{\phi}_{n,1}| < 10^{-5}$ for $m > 18$ which suggests that $M = 18$ is sufficient for generating the results with good accuracy in this case. However in the next section the profile for $\Omega_D$ is also included, which changes the

<figure><img src="image_9.png" /><figcaption>FIG. 9: Coefficients $\log(|\hat{\phi}_{n,m}|)$ of the function $C_m$ vs. $m$ for $M = 50$ points along $\psi$ in a computation for the mode $n = 15$ in case of temperature profile $T_4$ with $\psi_1 = 0.5$. A function $\log[A_0 \exp(-\beta m)]$ is plotted for $\beta = 0.7$ (red dashed line). $A_0$ is a constant.</figcaption></figure>

coefficients $\hat{\phi}_{n,m}$ values. Therefore to confirm the good accuracy for all other cases we consider $M = 50$ throughout this manuscript. Hereafter we will consider only the linear temperature profile $T_1$ and vary other important parameters for TIM instability.

## B. Variation in precession frequency

According to Eqs. (2,3), the precession frequency $\Omega_D$ depends on both $\psi$ and $\kappa$. In the previous section we consider a constant value of $\Omega_D = 1$. In this section we will first consider the effect of trapping parameter $\kappa$ on the precession frequency, and second, consider the $\psi$ dependency of $\Omega_D$ and investigate their effects on the TIM instability in the limit of $H_{eq} \approx E$.

<table><caption>TABLE III: Main input parameters used for studies in sec. IV B</caption><thead><tr><th>δ<sub>b,i</sub></th><th>ρ<sub>c,i</sub></th><th>C<sub>ad</sub></th><th>C<sub>pol</sub></th><th>T(ψ)</th><th>T<sub>eq</sub></th><th>n<sub>eq</sub></th><th>M</th><th>G<sub>T</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.03</td><td>0.1</td><td>0.1</td><td>T<sub>1</sub></td><td>1</td><td>1</td><td>50</td><td>0.25</td></tr></tbody></table>

### 1. $\kappa$ dependency of $\Omega_D$

In a banana orbit the critical poloidal angle $\theta_{crit}$, where $v_\parallel = 0$, is linked with the trapping parameter $\kappa$ as $\kappa^2 = \sin^2(\theta_{crit}/2)$. Therefore the acceptable values of $\kappa$ is $0 \le \kappa \le 1$. $\kappa = 0$ is associated with the particles having almost zero parallel velocity and therefore their motions are restricted close to the center of the banana and are called deeply trapped particles. $\kappa \sim 1$ is associated with the particles having maximum parallel velocity for
<figure><img src="image_10.png" /><figcaption>FIG. 10: Precession frequency $\bar{\omega}_d(\kappa)$ Eq. (3) for magnetic-shear $s_0 = 0.8$.</figcaption></figure>

the trapped particles, $\frac{|v_\parallel|}{|v_\perp|} \sim \sqrt{\frac{B_{max}}{B_{min}} - 1}$, whose reflecting point is situated near the position of maximum magnetic field $B_{max}$ at innerside of tokamak, and are called barely trapped particles. Eq. (3) presents the theoretical $\kappa$ dependency of the precession frequency $\bar{\omega}_d(\kappa)$. Fig. 10 presents this variation of $\bar{\omega}_d$ with $\kappa$ for a constant value of magnetic-shear $s_0 = 0.8$. In absence of $\psi$ dependency of precession frequency, we can write $\Omega_D(\kappa) = \bar{\omega}_d$ from Eq. (2). For a realistic safety factor profile [37], the magnetic shear $s_0(\psi)$ decreases almost linearly in $\psi$, and at $\psi = 0.5$ it has value $s_0 \sim 0.8$. In previous section-(IV A), we considered deeply trapped particles with $\kappa = 0$. Here we consider two different cases with same temperature profile $T_1$. In one case we consider constant value of precession frequency $\Omega_D = 0.6$ associated with barely trapped particles ($\kappa \sim 1$), and to compare the results we reconsider the previous case with $\Omega_D = 1$ ($\kappa = 0$) and temperature profile $T_1$. In second case to incorporate the effect of pitch-angle dependency in the TIM-instability we consider the entire $\Omega_D(\kappa)$ profile for the trapped particles. In this 2nd case the dispersion relation is calculated by doing the integration along $\kappa$ systemically according to Eq. (12), which was neglected for constant $\Omega_D$ by replacing $\int_0^1 \kappa K(\kappa^2) d\kappa = 1$. One important thing is that, for both cases we assume that the fraction of trapped particle $f_t \sim \sqrt{2\epsilon}$ remains constant, which is controlled by the parameter $C_{ad}$ in Eq. (12). If we consider radial dependency of inverse of aspect ratio as $\epsilon = r/R_0$, the fraction of trapped particle can be written as $f_t \sim (1-\psi)^{1/4}$. This justifies our approximation $f_t \sim \text{const.}$ throughout the region $\psi \in [0, 1]$. Therefore in 1st case with constant $\Omega_D = 0.6$, all trapped-particles are trapped near the separatrix (barely trapped), and for $\Omega_D = 1$, equal amount of particles are trapped near the center of battery, whereas in the 2nd case with entire $\Omega_D(\kappa)$ profile, equal amount of trapped particles are distributed over the entire trapped domain. The other essential parameters are taken from Tab. III.

<figure><img src="image_11.png" /><figcaption>FIG. 11: Growth rate $\gamma$ Vs. mode-number $n$ with $J_0$ from Eq. (4) and temperature $T_1(\psi)$, for constant values of precession frequency $\Omega_D = 1$ (solid blue line) and 0.6 (solid red line). The black solid line is for entire profile of $\Omega_D(\kappa)$ Eq. (3). Dots are the results from the linear TERESA simulations.</figcaption></figure>

Fig. 11 presents the growth-rates $\gamma$ for different values of mode numbers $n$ of TIM-instability for these cases. Since the TIM instability occurs due to the resonance of precession frequency $\Omega_D$ of particles with the wave frequency $\omega$, for low precession frequency $\Omega_D = 0.6$ the wave with smaller phase velocity (smaller $\omega$) resonates. In the equilibrium distribution $F_{eq}(E)$ Eq. (A1), there are large number of particles near lower velocity ($E \equiv \frac{1}{2}mv_{G\parallel}^2 + \mu B_G$) compared to higher velocity, as a consequence the charge separation due to $\nabla B$ drift in presence of temperature gradient would be higher near the lower velocity compared to higher velocity (higher $\Omega_D$), and generates stronger electric field which helps to enhance the density perturbation of the wave. Therefore $\Omega_D = 0.6$ has higher growth-rate compared to the case with higher precession frequency $\Omega_D = 1$. Since the real part of the TIM frequency $\omega_r \propto \Omega_D$ Eq. (A5), for lower $\Omega_D = 0.6$ value $\omega_r$ is significantly smaller compared to $\Omega_D = 1$ case. If we consider the entire profile of $\Omega_D(\kappa)$ for the constant magnetic-shear $s_0 = 0.8$ (Fig. 10), up to $\kappa \le 0.85$, the precession frequency $\Omega_D \ge 1$, and within $0.85 \le \kappa \le 1$ it has smaller value $\Omega_D < 1$. Since the equilibrium distribution $F_{eq}$ Eq. (A1) is independent of $\kappa$, there are equal number of particles at all $\kappa$ values. Therefore in this case, most of the particles (almost 80%) have precession frequency $\Omega_D > 1$, only 20% particles have precession frequency $\Omega_D < 1$. As a result the contribution from the $\Omega_D > 1$ dominates, therefore the waves with higher phase velocity (higher $\omega$) resonates with the particles, where less number of particles are available due to the Maxwellian particle distribution in $v_\parallel$. This decreases the growth-rate $\gamma$ compared to the $\Omega_D = 1$ case. The entire $\kappa$ dependent $\Omega_D(\kappa)$ profile decreases the growth-rate by 30–50% compared to constant $\Omega_D$ value at a constant $\kappa$ location. Dots are the growth-rates for different mode numbers $n$ from linear TERESA simulation.
<figure><img src="image_12.png" /><figcaption>FIG. 12: Precession frequency profiles $\Omega_D(\psi) = \Omega_{D0} + G_\Omega\psi$, with $G_\Omega = 0$ (solid red line), $G_\Omega = 0.25$ (solid blue line), $G_\Omega = 0.4$ (solid black line) and $G_\Omega = 0.6$ (solid magenta line). $\Omega_{D0} = 1$ for all the cases. The theoretical $\Omega_D(\psi)$ profile (solid cyan line) is calculated from Eq. (2) for $\kappa = 0$ and $q(r) = 1.1 + 2r^2$.</figcaption></figure>

## 2. $\psi$ dependency of $\Omega_D$

The theoretical dependency of precession frequency on radius $r$ is given by Eqs. (2-3) where both the safety factor $q(r)$, and magnetic shear $s_0(r) = \frac{r}{q(r)} \frac{dq}{dr}$ depends on $r$. For a particular $q(r)$ profile $q(r) = 1.1 + 2r^2$ profile [37], the magnetic flux-function is calculated from the integration $\psi(r) \sim -B_{min} \int_{r_{in}}^{a} \frac{r}{q(r)} dr$, and normalized according to Tab. I, where $L_\psi = |\psi(r_{in}) - \psi(a)|$ is the length of the simulation box in $\psi$ unit. Finally $r$ depends on $\psi$ as,

$$r = \sqrt{\frac{q(a)}{2} \left( \frac{q(r_{in})}{q(a)} \right)^{\psi} - \frac{q(0)}{2}}, \quad (31)$$

where $q(a)$, $q(0)$ and $q(r_{in})$ are the value of safety factor at $r = a$, $r = 0$ and $r = r_{in}$, where $a$ is the minor radius of tokamak and $r_{in}$ is the lower limit of $r$ integration in $\psi(r)$ expression, which helps to remove the singular nature of $\omega_d$ in Eq. (2) at $r \to 0$. Using this expression of $r$ Eq. (31) into $q(r)$ and $s_0(r)$, the $\psi$ dependent safety factor $q(\psi)$ and magnetic shear $s_0(\psi)$ can be obtained, and finally $\Omega_D(\psi, \kappa)$ can be calculated from Eqs. (2-3). The theoretical $\Omega_D(\psi)$ profile for $r_{in} = 0.2$ and $\kappa = 0$ is presented by solid cyan line in Fig. 12, which has almost constant value $\Omega_D \sim 1$ within $0 \le \psi \le 0.75$ and after that near the core region of tokamak ($\psi \to 1$), $\Omega_D$ increases abruptly with increase in $\psi$. For simplicity, first we consider $\Omega_D$ varies linearly with $\psi$ as:

$$\Omega_D(\psi) = \Omega_0 + G_\Omega \psi, \quad (32)$$

and study the effect of $\psi$ dependent $\Omega_D(\psi)$ on TIM instability. Here $\Omega_0 = 1$ is the value of $\Omega_D$ at $\psi = 0$, and $G_\Omega$ is the gradient in the $\Omega_D(\psi)$ profile. This linearly increasing $\Omega_D(\psi)$ profile for different values of $G_\Omega = 0$,

<figure><img src="image_13.png" /><figcaption>FIG. 13: Growth rate $\gamma$ Vs. mode-number $n$ with $J_0$ from Eq. (4) and temperature $T_1(\psi)$, and $\psi$ dependent precession frequency $\Omega_D(\psi) = \Omega_0 + 0.25\psi$ is presented in solid red line. Solid blue line presents the $\Omega_D = 1$ case. Dots are the results from the linear TERESA simulations.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 14: Potential solutions $|\phi_n(\psi)|/\phi_{max}$ for mode $n = 1$ with $J_0$ Eq. (4) for the precession frequency $\Omega_D = 1$ (dashed cyan), $\Omega_D = 0.6$ (solid red), $\Omega_D = 1 + 0.25\psi$ (dashed black), $\Omega_D(\kappa)$ (dashed blue).</figcaption></figure>

0.25, 0.4 and 0.6 are presented in Fig. 12. In this section the temperature profile $T_1(\psi)$ with $G_T = 0.25$, and $\Omega_D(\psi)$ profile with $G_\Omega = 0.25$ and $G_\Omega = 0$ Eq. (32) are taken into account. The other $\Omega_D(\psi)$ profiles will be considered in the next section Sec. V. The other essential parameters are taken from the Tab. III.

Fig. 13 presents the growth-rates of TIM modes $n$, for the $\Omega_D(\psi)$ profile Eq. (32) with $G_\Omega = 0.25$, and $G_\Omega = 0$. As discussed in the previous sub-section (IV B 1), for larger $\Omega_D$, the number of resonant particles decreases compared to the case with smaller $\Omega_D$ value, and as a consequence generates smaller growth-rate of the TIM instability. Since in the $\Omega_D$ profile with $G_\Omega = 0.25$ Eq. (32), for all $\psi$ value $\Omega_D \ge 1$ and the temperature profile $T_1(\psi)$ remains unchanged, the growth-rate $\gamma$ for all $n$ values is smaller compared to the case with $\Omega_D = 1$. Here
the dots present the results from the linear TERESA simulation. The real part of the frequency $\omega_r$ for the TIM instability is proportional to $\Omega_D$ Eq. (A5). Therefore in this case $\omega_r$ will be larger compared to the case with $\Omega_D = 1$. However we have not presented the $\omega_r - n$ profile within this manuscript. Fig. 14 presents the potential $\phi_n(\psi)$ profiles for the four different $\Omega_D$ profiles with constant temperature profile $T_1$, which are almost similar. Therefore in the limit $H_{eq} \approx E$ the potential profiles are almost independent of the precession frequency profiles. One important point is that, in case of $\psi$ dependent precession frequency $\Omega_D(\psi)$, the equilibrium Hamiltonian $H_{eq}$ depends on $\psi$ in more complicated manner, which was linear in $\psi$ for constant $\Omega_D$. Therefore the limit $H_{eq} \approx E$ is no more a good approximation, one has to consider the expression of exact Hamiltonian. We will discuss this issue in the next section.

V. EFFECT OF INVERSE GRADIENT LENGTH OF EQUILIBRIUM HAMILTONIAN $\kappa_\Lambda$ ON TIM INSTABILITY

In this section we consider the dispersion relation (Eq. (13-15)), which is derived from the exact Hamiltonian expression Eq. (7) in sec. II. In this case, according to Eq. (20), the threshold value of $\kappa_T$ for TIM instability, is higher compared to the case in the limit $H_{eq} \approx E$. Therefore in this section we have considered the temperature profile $T_1(\psi)$ with gradient $G_T = 3$ and 2 in Eq. (30) which was 0.25 in the previous cases Sec. IV in the limit $H_{eq} \approx E$. For investigating the effect of $\kappa_\Lambda$ on TIM instability, we consider precession frequency $\Omega_D(\psi)$ according to the theoretical expression Eq. (2-3) and also the simplified expression in the form of Eq. (32) with three different values of gradient $G_\Omega = 0, 0.4$ and $0.6$. Fig. 12 presents all these $\Omega_D(\psi)$ profiles. Considering

<table><caption>TABLE IV: Main input parameters used for studies in sec. V.</caption><thead><tr><th>δ<sub>b,i</sub></th><th>ρ<sub>c,i</sub></th><th>C<sub>ad</sub></th><th>C<sub>pol</sub></th><th>T(ψ)</th><th>T<sub>eq</sub></th><th>n<sub>eq</sub></th><th>M</th><th>G<sub>T</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.03</td><td>0.1</td><td>0.1</td><td>T<sub>1</sub></td><td>1</td><td>1</td><td>50</td><td>3, 2</td></tr></tbody></table>

the theoretical expression of $\Omega_D(r, \kappa)$ Eq. (2-3) and the safety factor profile as $q(r) = 1.1 + 2r^2$, the expression of $\Lambda_D(r, \kappa)$ can be written as,

$$ \Lambda_D(r, \kappa) = 1 + \frac{4}{q(a) \ln\left(\frac{q(a)}{q(r_{in})}\right)} \left[ -\lambda_2(\kappa)r + 4\sqrt{2q(0)} \right. \\ \left. \lambda_1(\kappa) \tan^{-1}\left(\sqrt{\frac{2}{q(0)}}r\right) + \lambda_2(\kappa)\sqrt{\frac{1}{2}(q(a) - q(0))} \right] \quad (33) \\ - 4\sqrt{2q(0)}\lambda_1(\kappa) \tan^{-1}\left(\sqrt{\frac{q(a)}{q(0)}} - 1\right) $$

<figure><img src="image_15.png" /><figcaption>FIG. 15: Different profiles of $\Lambda_D(\psi) = 1 + \int \Omega_D(\psi)d\psi$, associated with three different $\Omega_D(\psi)$ profiles. Solid red, black and magenta lines are associated with $G_\Omega = 0, 0.4$ and $0.6$, respectively in Eq. (32). The solid cyan line presents $\Lambda_D(\psi)$ profile for the theoretical expression of $\Omega_D(\psi)$ Eq. (2-3).</figcaption></figure>

where $\lambda_1(\kappa) = \left(\frac{\varepsilon(\kappa^2)}{K(\kappa^2)} + \kappa^2 - 1\right)$, and $\lambda_2(\kappa) = \left(10\frac{\varepsilon(\kappa^2)}{K(\kappa^2)} + 8\kappa^2 - 9\right)$. Substituting the value of $r$ from Eq. (31), $\psi$ dependent $\Lambda_D(\psi, \kappa)$ can be calculated. This $\Lambda_D(\psi)$ profile Eq. (33) with $\kappa = 0$ is presented in Fig. 15 with solid cyan line. The other profiles of $\Lambda_D(\psi)$ in Fig. 15 are, for the simplified $\Omega_D(\psi)$ profile in Eq. (32) with $G_\Omega = 0, 0.4$ and $0.6$. Therefore the $\Lambda_D(\psi)$ profile for the simplified $\Omega_D(\psi)$ expression Eq. (32) with $G_\Omega = 0$, is almost similar with the $\Lambda_D(\psi)$ profile for the theoretical expression of $\Omega_D(\psi, \kappa)$ with $\kappa = 0$. The values of $\kappa_\Lambda = \frac{\Omega_D}{\Lambda_D}$ for all the four profiles at $\psi = 0$ is $\kappa_\Lambda = 1$, and then decreases with different rates as $\psi$ increases. Here first, we consider the temperature profile $T_1$ with gradient $G_T = 3$ and study the TIM instability for all the four $\Lambda_D$ and $\Omega_D$ profiles. After that, for understanding the effect of temperature gradient on TIM instability in this new modified model with exact equilibrium Hamiltonian, we decrease the temperature gradient $G_T = 2$ and consider the theoretical expression of $\Omega_D$ as Eqs. (2-3) and $\Lambda_D$ profile as Eq. (33) for $\kappa = 0$. Other essential parameters are taken from the tab. IV. The potential solution $\phi_n(\psi)$ is obtained by solving the differential equation Eq. (16) with the expression for $N_n$ Eq. (13), using the spectral method. Due to high temperature gradient $G_T = 3$ and $2$, the value of $T(\psi)$ is very high compared to previous cases, which makes the particle equilibrium distribution $F_{eq}(\psi, \kappa, E)$ Eq. (10) broaden to high $E$ value. Therefore in this case we have to increase the maximum limit of $E$ as $E \in [0, 45]$ for the numerical integration along $E$ direction, whereas for the previous cases with $G_T = 0.25$ we considered $E \in [0, 20]$.

Fig. 16 presents the real part of the TIM instability for different mode numbers $n$ for all these five cases. According to Eq. (18), the case with higher $\kappa_\Lambda$ value has higher $\omega_r$, therefore the case with larger $G_\Omega$ has higher
<figure><img src="image_16.png" /><figcaption>FIG. 16: Real part of the frequency $\omega_r$ of TIM instability for $T_1$ with $G_T = 3$ and $\Omega_D(\psi)$ profiles according to the theoretical expression Eq. (2-3) with $\kappa = 0$ (solid cyan line), and from the simplified expression Eq. (32) with $G_\Omega = 0$ (solid red line), $G_\Omega = 0.4$ (solid blue line) and $G_\Omega = 0.6$ (solid black line). The solid magenta line: for $G_T = 2$ and the theoretical expression of $\Omega_D(\psi)$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 17: Growth rate $\gamma$ of different modes $n$, $T_1$ with $G_T = 3$ and different $\Omega_D(\psi)$ profiles, solid red, blue and black lines are for simplified $\Omega_D(\psi)$ expression Eq. (32) with $G_\Omega = 0$, $G_\Omega = 0.4$ and $G_\Omega = 0.6$ respectively. Solid cyan line: $\Omega_D(\psi)$ theoretical expression Eq. (2-3) with $\kappa = 0$, $G_T = 3$ and solid magenta line: for $G_T = 2$ and the theoretical expression of $\Omega_D(\psi)$. The dots are the results from linear TERESA simulations.</figcaption></figure>

$\omega_r$. In the theoretical expression of $\Omega_D$ for $\kappa = 0$ the value of $\Omega_D \sim 1$ for $\psi \in [0, 0.75]$ and only for $\psi \ge 0.75$, it has values greater than unity. Therefore the $\omega_r$ profile of TIM instability, in this case is very close to the case with $\Omega_D = 1$ profile. Keeping fixed the $\Omega_D$ profile, if we decrease the $\kappa_T$ value by decreasing $G_T = 2$, the $\omega_r$ value decreases with a very small amount. Fig. 17 presents the growth-rate $\gamma$ for different mode numbers $n$ of TIM instability. Since $\kappa_\Lambda$ in the expression of $\mathcal{N}$ Eq. (13) reduces the effect of $\kappa_T$, the growth rates in all these cases are very smaller compared to the previous

<figure><img src="image_18.png" /><figcaption>FIG. 18: Potential solutions $|\phi_n(\psi)|/\phi_{max}$ for mode $n = 1$, $T_1$ with $G_T = 3$ and different $\Omega_D(\psi)$ profiles. Solid red, blue and black lines are for simplified $\Omega_D(\psi)$ expression Eq. (32) with $G_\Omega = 0$, $G_\Omega = 0.4$ and $G_\Omega = 0.6$ respectively. Solid cyan line: $\Omega_D(\psi)$ theoretical expression Eq. (2-3) with $\kappa = 0$, Solid magenta line: for $G_T = 2$ and the theoretical expression of $\Omega_D(\psi)$.</figcaption></figure>

cases with $\kappa_\Lambda = 0$ and $\Lambda = 1$. Moreover the increase in $\kappa_\Lambda$ decreases the $\gamma$ value. Among the three $\Omega_D$ profiles, the cases with highest $G_\Omega = 0.6$ value has higher $\kappa_\Lambda$ value, which gives the lowest growth-rate. As in the theoretical expression of $\Omega_D$ for $\kappa = 0$ has value $\Omega_D \sim 1$ within the region $\psi \in [0, 0.75]$ and after that it increases, therefore the growth-rate $\gamma$ in this case is almost similar to the case with $\Omega_D = 1$ and have slightly smaller value due to $\Omega_D \ge 1$ within the region $\psi \ge 0.75$. Therefore we can conclude that the complicated theoretical expression of $\Omega_D$ for $\kappa = 0$ can be simplified as $\Omega_D = 1$ with a good accuracy for studying the TIM instability. If we decrease the temperature gradient $G_T = 2$, for the same theoretical $\Omega_D$ profile, the growth rate $\gamma$ decreases by a significant amount ($\sim 70\%$) due to decrease in $\kappa_T$ value. The dots are the results from the linear TERESA simulations for all cases. Fig. 18 presents the potential profiles $\phi_n(\psi)$ for $n = 1$ and for all the five different cases with different $\Omega_D(\psi)$ profiles. The $\phi_n(\psi)$ solution are different for different $\Omega_D(\psi)$ profiles, whereas the $\phi_n(\psi)$ profiles for different $\Omega_D$ profiles were almost similar for the limit $H_{eq} \approx E$. In this new model $\phi_n(\psi)$ profiles are varying with both the temperature profiles and the precession frequency profiles.

Throughout this paper in the linear analysis, the non-linear term $[J_{0,s}\dot{\phi}, \hat{f}_s]_{\alpha,\psi}$ is neglected, by considering that $\hat{f}_s$ is very small and the equilibrium distribution function $F_{eq,s}$ remains unchanged during the time evolution. However we have verified that, these assumptions are valid only for small amplitude potential $|e\phi/T| \ll 1$, and as the amplitude of the potential grows beyond a certain value ($e\phi/T \sim 0.1$), due to strong nonlinear wave-particle and/or wave-wave interactions, the growth rate $\gamma(t)$ decreases and finally nonlinear saturation occurs. Since the
nonlinear evolution of the TIM is not the focus of this
study, we have not presented those results here.

VI. CONCLUSIONS

In this work we proposed an alternative way to solve
reduced gyro-bounce averaged kinetic model for trapped
particle dynamics within linear limit using a spectral
method [28, 29]. This method is computationally very
fast compared to the semi-Lagrangian method based
solver TERESA [8, 25, 26]. Using this method we have
investigated the trapped ion mode (TIM) instability.
Unlike the local linear analysis [9] of trapped particle
mode instability, our proposed method can incorporate
the entire profiles (0 ≤ $\psi$ ≤ 1) of all the essential pa-
rameters, and is not restricted by their local values at
$\psi$ = 0. In that respect our method is a global-linear
analysis of trapped particle modes instability. Also the
dependency of trapped particle drift velocity on mag-
netic poloidal flux function $\psi$, is newly incorporated in
the gyro-bounce averaged trapped particle model by con-
sidering the exact expression of equilibrium Hamilto-
nian $H_{eq}(\psi, \kappa) = E\Lambda_D(\psi, \kappa)$ in quasi-neutrality equation
and gyro-bounce averaged equilibrium distribution func-
tion $F_{eq}$, which were previously simplified in the limit
$H_{eq} \approx E$. With this new Hamiltonian, a new quantity
$\kappa_\Lambda$, that measures the inverse gradient length of equilib-
rium Hamiltonian, appears in the dispersion relation of
the TIM instability. All the previous results in the limit
$H_{eq} \approx E$ can be recovered by substituting $\Lambda_D = 1$ and
$\kappa_\Lambda = 0$. The quantity $\kappa_\Lambda$ reduces the effect of $\kappa_T$ and
$\kappa_n$, and as a consequence reduces the growth-rate of TIM
instability.

In tokamak plasma, trapped particle modes (resonant-
branch) are driven by the resonant interaction with the
precession motion of trapped particles, and these modes
become unstable in presence of density inhomogeneity,
gradient in magnetic field, and above a critical gradient
of temperature. Therefore the effect of different tem-
perature profiles T(ψ) and precession-frequency Ω<sub>D</sub> on
the linear TIM instability are investigated. The results
for all the cases are compared with the linear TERESA
simulations. First we consider the model in the limit
H<sub>eq</sub> ≈ E and validate our spectral method based solver
with the TERESA simulation for the expression of gyro-
bounce average operator J<sub>0</sub> = 1. The solutions from this
J<sub>0</sub> = 1 case is used during calculation of gyro-bounce
average of potential J<sub>0</sub>φ using Padé expression, for the
mode number n = 1 at the 1<sup>st</sup> iteration. For studying
the effect of temperature gradient on TIM instability, we
vary the temperature profiles which are relevant to dif-
ferent regions of a tokamak plasma experiment and keep
fixed the normalized precession frequency value Ω<sub>D</sub> = 1.
Depending on the inverse temperature gradient length
κ<sub>T</sub> value, the growth-rate of TIM instability for different
temperature profiles is different. The profile with higher
κ<sub>T</sub> value gives higher growth rate, and the highest grow-

ing mode also shifted towards higher mode number. The
real part of frequency $\omega_r - n$ profile for all the tempera-
ture profiles with a constant $\Omega_D$ value, follow a relation
$\omega_r \sim 2n\Omega_D$. Then we have studied the effect of variations
in precession frequency profile $\Omega_D$ on TIM instability in
the limit $H_{eq} \approx E$. For lower value of $\Omega_D$ the waves
with smaller velocity make resonance with the trapped
ion motion. Due to Maxwellian energy distribution of
particles, a large number of particles near lower veloc-
ity takes parts in the instability generation mechanism,
and yields stronger instability for smaller $\Omega_D$. Hence
we get larger growth rate for smaller $\Omega_D$. The effect of
pitch-angle dependency of trapped particles on the TIM
instability is investigated by considering the entire $\kappa$ de-
pendent profile of $\Omega_D$. Since for a constant magnetic
shear $s_0 = 0.8$, for most of the $\kappa$ value $\Omega_D(\kappa) > 1$ except
near the separatrix $0.85 \le \kappa \le 1$, the growth-rate of TIM
instability is smaller compared to $\Omega_D = 1$ case. The $\psi$
dependency of the precession frequency $\Omega_D(\psi)$ on TIM
instability in the limit $H_{eq} \approx E$, is also investigated for
a simplified linearly increasing $\Omega_D$ profile. Since for this
profile the values of precession frequency $\Omega_D > 1$ for all
the $\psi > 0$ location, the growth-rate of the TIM instabil-
ity is smaller than the case with $\Omega_D = 1$. In the limit
$H_{eq} \approx E$, the potential profiles are almost independent
of different $\Omega_D(\psi, \kappa)$ profiles.

The effect of $\psi$ dependency of $\Omega_D$ on the TIM insta-
bility in the newly modified model with exact equilib-
rium Hamiltonian expression is investigated by consider-
ing a theoretical expression of $\Omega_D$ and a simplified lin-
early increasing function of $\Omega_D$. For $\kappa = 0$ the theo-
retical expression of $\Omega_D$ has value $\Omega_D \sim 1$ within the
region $0 \le \psi \le 0.75$, and after that it increases. The
growth-rate $\gamma$ of TIM instability for the particular theo-
retical $\Omega_D$ profile is almost similar to the $\Omega_D = 1$ case,
and due to $\Omega_D \ge 1$ with in the region $0.75 \le \psi \le 1$,
it is slightly smaller than $\Omega_D = 1$ case. With increase
of the slope $G_\Omega$ in the simplified linearly increasing $\Omega_D$
profile, which increases $\kappa_\Lambda$ value, the growth-rate of dif-
ferent modes $n$ for the instability decreases. Moreover
after the new modification, the decrease in slope $G_T$ of
the temperature profile (decrease $\kappa_T$ value) for a fixed
$\Omega_D$ profile, decreases the growth rate, which is consis-
tent with the results in the limit $H_{eq} \approx E$. In this newly
modified model the potential solutions are different for
different $\Omega_D$ profiles.

Acknowledgements

This work was funded by the Agence Nationale de
la Recherche for the Project GRANUL (ANR-19-CE30-
0005). This work was granted access to the HPC re-
sources of EXPLOR (Project No. 2017M4XXX0251),
and of CINECA MARCONI under Project FUA35
GSNTITE. This work has been carried out within the
framework of the EUROfusion Consortium, funded by
the European Union via the Euratom Research and
Training Programme (Grant Agreement No 101052200-EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. We are grateful to the anonymous reviewers for their constructive remarks.

## Appendix A: Linear analysis within the limit

$$H_{eq} \approx E$$

The global linear analysis for the case with the limit $H_{eq} \approx E$ can be derived in the similar way, as discussed in Sec. II 1 Here the normalized equilibrium distribution $F_{eq,s}$ is independent of $(\alpha, \kappa, t)$, and has the form of a two dimensional Maxwellian energy distribution function

$$F_{eq,s}(\psi, E) = \frac{n_s(\psi)}{T_s^{3/2}(\psi)} \exp\left(-\frac{E}{T_s(\psi)}\right). \quad (A1)$$

After substituting $F_{eq,s}$ from Eq. (A1), $\tilde{f}_s = \sum_{n,\omega} f_{s,n,\omega}(\psi, E, \kappa) \exp\{i(n\alpha - \omega t)\}$ and $\tilde{\phi} = \sum_{n,\omega} \phi_{n,\omega}(\psi) \exp\{i(n\alpha - \omega t)\}$, in Eq. (9), the solution of Vlasov equation in Fourier space become,

$$f_{s,n,\omega}(\psi, E, \kappa) = \frac{n\left[\kappa_n(\psi) + \kappa_T(\psi)\left(\frac{E}{T_s(\psi)} - \frac{3}{2}\right)\right]}{Z_s^{-1}n\Omega_D(\kappa)E - \omega} \quad (A2) \\ \{J_{0,n,s}\phi_{n,\omega}(\psi)\} F_{eq,s}(\psi, E),$$

In this case the elementary volume in phase-space can be written as $d^3v = 4\pi\sqrt{2}m^{-3/2}\sqrt{E}dE \frac{d\lambda}{4\Omega_D}$. Using this volume element $d^3v$ and $f_{n,\omega}$ from Eq. (A2) the expression of $N_s$ in Eq. (12) can be written (in the limit of a constant pitch-angle) as,

$$N_{n,s}(\psi) = \frac{1}{n_{eq}} \int_{0}^{\infty} J_{0,n,s} \left[ \frac{\kappa_n(\psi) + \kappa_T(\psi) \left( \frac{E}{T_s(\psi)} - \frac{3}{2} \right)}{Z_s^{-1} \Omega_D (E - \chi_s)} \right. \quad (A3) \\ \left. \{J_{0,n,s}\phi_{n,\omega}(\psi)\} \frac{n_s(\psi)}{T_s^{3/2}(\psi)} \frac{\exp\left(-\frac{E}{T_s(\psi)}\right)}{\sqrt{E}dE}, \right]$$

where $\chi_s = \frac{\omega}{nZ_s^{-1}\Omega_D}$. The expression of $C_n$ Eq. (14) remains unaltered. Therefore the dispersion relation becomes,

$$C_n \phi_{n,\omega} = N_{n,i}^* \phi_{n,\omega} - N_{n,e}^* \phi_{n,\omega} \quad (A4)$$

where $N_{n,s}^* = \frac{N_{n,s}}{\phi_{n,\omega}}$. The final form of the 2nd-order differential equation for $\phi_n(\psi)$ Eq. (16) remains unchanged, except $N_n^*$ inside $Q_n(\psi)$ Eq. (17), where $N$ is given by Eq. (A3). Comparing the expressions of $N$ from Eq. (A3) and Eq. (13), the main difference in the dispersion relation arises due to the absence of the term $-\frac{E\Lambda_D}{T}\kappa_\Lambda$ in $N$ in the limit $H_{eq} \approx E$, which reduces the effective contribution from $\kappa_n$ and $\kappa_T$ in case of exact Hamiltonian. Therefore the growth-rate of the TIM instability in the limit $H_{eq} \approx E$, is significantly higher compared to the case with exact Hamiltonian.

As discussed in Sec. II 2 the local linear stability analysis of this reduced gyro-bounce averaged model in the limit $H_{eq} \approx E$ can be performed and the threshold frequency value of TIM instability can be derived as [9]:

$$\omega_r = \frac{(\frac{3}{2}\kappa_T - \kappa_n)}{\kappa_T} n\Omega_D T_0, \quad (A5)$$

where $T_0$ is the temperature at $\psi = 0$. The threshold value of $\kappa_T$ for the TIM instability can be written as

$$\kappa_{T,th} = \frac{C_n \Omega_{D0}}{\int_{0}^{\infty} J_{0,n}^2 \exp(-\xi) \sqrt{\xi} d\xi} \quad (A6)$$

where $\xi = \frac{E}{T}$.

[1] A. M. Dimits, B. I. Cohen, W. M. Nevins, and D. E. Shumaker, Nucl. Fusion 41, 1725 (2001).

[2] P. N. Guzdar, L. Chen, W. M. Tang, and P. Rutherford, Phys. Fluids 26, 3 (1983).

[3] M. Ottaviani, W. Horton, and M. Erba, Plasma Phys. Control. Fusion 39, 1461 (1997).

[4] B. W. Stallard, C. M. Greenfield, G. M. Staebler, C. L. Retting, and et. al., Phys. Plasmas 6, 1978 (1999).

[5] T. Dannert and F. Jenko, Phys. Plasmas 12, 072309 (2005).

[6] D. L. Brower, W. A. Peebles, S. K. Kim, N. C. Luhmann, W. M. Tang, and P. E. Phillips, Phys. Review Let 59, 48 (1987).

[7] H. Eubank, R. Goldston, V. Arunasalam, M. Bitter, and et. al., Phys. Rev. Letter 43, 270 (1979).

[8] G. Depret, X. Garbet, P. Bertrand, and A. Ghizzo, Plasma Phys. Control. Fusion 42, 949 (2000).

[9] T. Drouot, E. Gravier, T. Reveille, A. Ghizzo, P. Bertrand, X. Garbet, Y. Sarazin, and T. Cartier-Michaud, Eur. Phys. J. D 68, 280 (2014).

[10] X. Garbet, Plasma Phys. Control. Fusion 43, A251 (2001).

[11] A. Ghizzo and F. Palermo, Phys. plasmas 22, 082304 (2015).

[12] H. Chen and L. Chen, Phys. Rev. Let. 128, 025003 (2022).

[13] T. Drouot, E. Gravier, T. Reveille, M. Sarrat, M. Collard, and et. al., Phys. Plasmas 22, 082302 (2015).
[14] R. E. Waltz, G. D. Kerbel, and J. Milovich, Phys. Plasmas **1**, 7 (1994).

[15] M. Ottaviani, F. Romanelli, R. Benzi, M. Briscolini, P. Santangelo, and S. Succi, Phys. Fluids **B** **2**, 1 (1990).

[16] W. Dorland and G. W. Hammett, Phys. Fluids **5**, 812 (1993).

[17] Z. Lin and L. Chen, Phys. Plasmas **8**, 1447 (2001).

[18] R. D. Sydora, V. K. Decyk, and J. M. Dawson, Plasma Phys. Control. Fusion **38**, A281 (1996).

[19] W. W. Lee, Journal Computational Phys. **72**, 243 (1987).

[20] Y. Idomura, S. Tokuda, and Y. Kishimoto, Nucl. Fusion **43**, 234 (2003).

[21] G. Manfredi, M. Shoucri, R. O. Dendy, A. Ghizzo, and P. Bertrand, Phys. Plasmas **3**, 202 (1996).

[22] J.-M. Kwon, L. Qi, S. Yi, and T.S.Hahmba, Computer Phys. Commun. **215**, 81 (2017).

[23] B. H. Fong and T. S. Hahm, Phys. Plasmas **6**, 188 (1999).

[24] A. Ghizzo, M. E. Mouden, D. D. Sarto, X. Garbet, and Y. Sarazin, Transport Theory Statistical Phys. **40**, 382 (2011).

[25] Y. Sarazin, V. Grandgirard, E. Fleurence, X. Garbet, P. Ghendrih, P. Bertrand, and G. Depret, Plasma Phys. Control. Fusion **47**, 1817 (2005).

[26] G. Darmet, P. Ghendrih, Y. Sarazin, X. Garbet, and V. Grandgirard, Commun. Nonlinear Sc. and Numerical

Simulation **13**, 53 (2008).

[27] M. Lesur, C. Djerroud, K. Lim, E. Gravier, M. Idouakass, and et. al., Nucl. Fusion **60**, 036016 (2020).

[28] E. Gravier, E. Plaut, X. Caron, and M. Jenny, Eur. Phys. J. D **67**, 7 (2013).

[29] E. Gravier and E. Plaut, Phys. Plasmas **20**, 042105 (2013).

[30] T. Drouot, Ph.D. thesis, Université de Lorraine (2015).

[31] J. Plemelj, Monatshefte für Mathematik und Physik **19**, 205 (1908).

[32] V. G. Pryimak and T. Miyazaki, J. Comput. Phys **142**, 370 (1998).

[33] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, SIAM J. Opt. **9**, 112 (1998).

[34] M. Lesur, T. Cartier-Michaud, T. Drouot, P. H. Dia mond, and et. al., Phys. Plasmas **24**, 012511 (2017).

[35] J. Jacquinot and the JET team, Plasma Phys. Control. Fusion **41**, A13 (1999).

[36] R. Pánek, J. Adámek, M. Aftanas, P. Bilková, and et. al., Plasma Phys. Control. Fusion **58**, 014015 (2015).

[37] M. Gobbin, D. Bonfiglio, A. H. Boozer, A. W. Cooper, D. F. Escande, S. P. Hirshman, and et. al., Phys. Plasmas **18**, 062505 (2011).