2.1. Linear analysis

We linearize the Vlasov and quasi-neutrality equations by splitting each quantity into an equilibrium and a perturbed component, such that

(2.3)\begin{gather} f_s=f_{0,s}+f_{1,s}, \end{gather}

(2.4)\begin{gather} \dot{\boldsymbol{R}}=\dot{\boldsymbol{R}}_0+\dot{\boldsymbol{R}}_1, \end{gather}

(2.5)\begin{gather} \dot{E}=\dot{E}_0+\dot{E}_1, \end{gather}

(2.6)\begin{gather} \varPhi=\varPhi_1, \end{gather}

where $\dot {\boldsymbol {R}}_0$ is the unperturbed particle velocity, i.e. $\dot {\boldsymbol {R}}_0= \boldsymbol {v}_{\|}+\boldsymbol {v}_{\boldsymbol {\nabla } B}+\boldsymbol {v}_{curv B}$, and $E_0=({m(v_{\|}^2+v_{\perp }^2)})/{2}$.

2.2. Linear Vlasov equation

Substituting (2.3)–(2.6) in (2.1), and neglecting second- and higher-order terms, the linear Vlasov equation reads

(2.7)\begin{equation} \frac{\partial f_{1,s}}{\partial t}+\dot{\boldsymbol{R}}_0\boldsymbol{\cdot}\frac{\partial f_{1,s}}{\partial\boldsymbol{R} }={-}\dot{E}_1\frac{\partial f_{0,s}}{\partial E }. \end{equation}

This linear Vlasov equation can be further simplified by splitting the perturbed distribution function into an adiabatic and a non-adiabatic component:

(2.8)\begin{equation} f_{1,s}=E_1\frac{\partial f_{0,s}}{\partial E}+h_s. \end{equation}

Substituting (2.8) in (2.7), we obtain the equation for the non-adiabatic part of the perturbed distribution function:

(2.9)\begin{equation} \frac{\partial h_s}{\partial t}+\dot{\boldsymbol{R}}_0\boldsymbol{\cdot}\frac{\partial h_s}{\partial\boldsymbol{R} }={-}\dot{E}_1\frac{\partial f_{0,s}}{\partial E }. \end{equation}

Using the expression of the equilibrium velocity and perturbations of the form $x \longrightarrow \exp {({\rm i}k_r-{\rm i}\omega t)}+{\rm c.c.}$ in (2.9), we obtain the expression

(2.10)\begin{equation} \left(\omega_{t,s}\frac{\partial}{\partial \theta}-{\rm i}(\omega+\omega_{d,s})\right)h_s={\rm i}\omega E_1\frac{\partial f_{0,s}}{\partial E}, \end{equation}

where $\omega _t={v_{\|}}/{qR_0}$ is the transit frequency, $q$ is the safety factor, $\omega _{d,s}=\bar {\omega }_{d,s}\sin {\theta }$ is the drift frequency, with $\bar {\omega }_{d,s}=({cm_sk_r}/{Z_seB_0R_0})(v_{\|}^2+{v_{\perp }^2}/{2})$, $E_1=Z_seJ_{0,s}\varPhi _1$ and $k_r$ is the radial wavenumber. Since GAMs are predominantly zonal, we can further divide the non-adiabatic part of the perturbed distribution function into a zonal and a non-zonal part:

(2.11)\begin{equation} h=\bar{h}+\delta h. \end{equation}

Similarly, we write the scalar potential as

(2.12)\begin{equation} \varPhi_1=\bar{\varPhi}+\tilde{\varPhi}, \end{equation}

where the overbar represents the zonal components. Using these definitions, and making a flux surface average of the gyro-kinetic equation to eliminate the zonal component of the non-adiabatic part of the perturbed distribution function, the linear Vlasov equation reduces to

(2.13)\begin{equation} \left(\omega_{t,s}\frac{\partial}{\partial \theta}-{\rm i}(\omega+\omega_{d,s})\right)\delta h_s={\rm i}Z_s\frac{\partial f_{0,s}}{\partial E}\left(\omega J_{0,s}\tilde{\varPhi}-\omega_{d,s} J_{0,s}\bar{\varPhi}\right). \end{equation}

The corresponding vorticity equation is obtained by multiplying this relation by the gyro-average operator and integrating over the velocity space:

(2.14)\begin{equation} \left\langle J_{0,s}\left(\omega_{t,s}\frac{\partial}{\partial \theta}-{\rm i}(\omega+\omega_{d,s})\right)\delta h_s\right\rangle_W=\left\langle {\rm i}Z_s\frac{\partial f_{0,s}}{\partial E}\left(\omega J_{0,s}^2\tilde{\varPhi}-\omega_{d,s} J_{0,s}^2\bar{\varPhi}\right)\right\rangle_W. \end{equation}

Considering all the changes of variable we have made, the perturbed distribution function has the form

(2.15)\begin{equation} f_{1,s}=Z_seJ_{0,s}\frac{\partial f_{0,s}}{\partial E}\tilde{\varPhi}+\delta h_s. \end{equation}

2.3. Linear quasi-neutrality equation

Using a similar approach, we substitute (2.3) and (2.6) into (2.2) and we thus obtain the following equation:

(2.16)\begin{equation} \dfrac{m_ic^2}{B^2}k_r^2\bar{\varPhi}=e^2\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W\left(1+\dfrac{\left\langle J_{0,e}^2\dfrac{\partial f_{0,e}}{\partial E}\right\rangle_W}{\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W}\right)\tilde{\varPhi}+e\left\langle J_{0,i}\delta h_i\right\rangle_W, \end{equation}

where $\left \langle \cdots \right \rangle _W$ represents the integral over velocity space. The non-adiabatic part of the perturbed electron distribution function has been neglected in accordance with our assumptions.

2.4. Ordering of the gyro-kinetic equation

The gyro-kinetic equation (2.13) describes a wide range of phenomena at different time scales. In order to study GAMs, we need to apply an appropriate ordering that will filter out time scales which are irrelevant to GAM dynamics. The GAM frequency is of the order of ion sound frequency. The ordering is done by comparing this frequency with the characteristic frequencies in our system, i.e. $\omega _{t,s}$, $\omega _{d,s}$:

(2.17a–c)\begin{gather} \dfrac{\omega_{t,i}}{\omega} \sim O(1), \quad \frac{\omega_{d,i}}{\omega}\sim O(\epsilon), \quad \frac{\delta h_i}{f_{0,i}} \sim O(\epsilon), \end{gather}

(2.18a–c)\begin{gather} \frac{\omega_{t,e}}{\omega} \gg 1, \quad \frac{\omega_{d,e}}{\omega} \sim O(1), \quad \frac{\delta h_e}{f_{0,e}} \sim O(\epsilon). \end{gather}

To leading order, the ion and electron gyro-kinetic equations are, respectively,

(2.19)\begin{gather} \left(\frac{\omega_{t,i}}{\omega}\frac{\partial}{\partial \theta}-{\rm i}\right)\delta h_i={\rm i}Z_ie\frac{\partial f_{0,i}}{\partial E}\left(J_{0,i}\tilde{\varPhi}-J_{0,i}\frac{\omega_{d,i}}{\omega}\bar{\varPhi}\right), \end{gather}

(2.20)\begin{gather} \delta h_e=0. \end{gather}

2.5. General form of dispersion relation

We consider the following form for the non-zonal perturbed ion distribution function and scalar potential:

(2.21)\begin{gather} \tilde{\varPhi}=\tilde{\varPhi}_s\sin{\theta}+\tilde{\varPhi}_c\cos{\theta}, \end{gather}

(2.22)\begin{gather} \delta h_i=\delta h_{i,s}\sin{\theta}+\delta h_{i,c}\cos{\theta}. \end{gather}

Substituting these relations into (2.19) and separating the sine and cosine components, we obtain

(2.23)\begin{gather} \delta h_{i,s}=\dfrac{{\rm i}Z_iJ_{0,i}\dfrac{\partial f_{0,i}}{\partial E}}{\left(\dfrac{\omega_{t,i}}{\omega}\right)^2-1} \left[-{\rm i}\left(\tilde{\varPhi}_s-\left(\dfrac{\bar{\omega}_{d,i}}{\omega}\right)\bar{\varPhi}\right)+\dfrac{\omega_{t,i}}{\omega}\tilde{\varPhi}_c\right], \end{gather}

(2.24)\begin{gather} \delta h_{i,c}={-}\dfrac{{\rm i}Z_iJ_{0,i}\dfrac{\partial f_{0,i}}{\partial E}}{\left(\dfrac{\omega_{t,i}}{\omega}\right)^2-1} \left[\dfrac{\omega_{t,i}}{\omega}\left(\tilde{\varPhi}_s -\left(\dfrac{\bar{\omega}_{d,i}}{\omega}\right)\bar{\varPhi}\right)+{\rm i}\tilde{\varPhi}_c\right]. \end{gather}

Following the same procedure with the quasi-neutrality equation:

(2.25)\begin{gather} \tilde{\varPhi}_c={-}\dfrac{\left\langle J_{0,i}\delta h_{i,c}\right\rangle_W}{e\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W\left(1+\dfrac{\left\langle J_{0,e}^2\dfrac{\partial f_{0,e}}{\partial E}\right\rangle_W}{\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W}\right)}, \end{gather}

(2.26)\begin{gather} \tilde{\varPhi}_s={-}\dfrac{\left\langle J_{0,i}\delta h_{i,s}\right\rangle_W}{e\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W\left(1+\dfrac{\left\langle J_{0,e}^2\dfrac{\partial f_{0,e}}{\partial E}\right\rangle_W}{\left\langle J_{0,i}^2\dfrac{\partial f_{0,i}}{\partial E}\right\rangle_W}\right)}. \end{gather}

Taking the flux surface average of the quasi-neutrality equation (2.19) and the vorticity equation (2.14), we have, respectively,

(2.27)\begin{gather} \frac{m_ic^2}{B^2}k_r^2\bar{\varPhi}=e\overline{\left\langle J_{0,i}\delta h_i\right\rangle}_W, \end{gather}

(2.28)\begin{gather} \overline{\left\langle J_{0,i}\delta h_i\right\rangle}_W={-}\frac{\overline{\left\langle J_{0,i}\omega_{d,i}\delta h_i\right\rangle}_W}{\omega}. \end{gather}

Substituting (2.28) into (2.27) and evaluating the flux surface average, we obtain

(2.29)\begin{equation} \frac{2m_i}{B^2}k_r^2\bar{\varPhi}={-}\frac{e}{c^2\omega}\left\langle J_{0,i}\bar{\omega}_{d,i}\delta h_{i,s}\right\rangle. \end{equation}

Substituting (2.23), we obtain the general dispersion relation of GAMs to leading order:

(2.30)\begin{equation} \dfrac{2m_ic^2k_r^2}{B^2}\bar{\varPhi}={-}e^2\omega\left\langle \dfrac{J_{0,i}^2\bar{\omega}_{d,i}\dfrac{\partial f_{0,i}}{\partial E}}{\omega_{t,i}^2-\omega^2}\right\rangle_W\tilde{\varPhi}_s+e^2\left\langle \dfrac{J_{0,i}^2\bar{\omega}_{d,i}^2\dfrac{\partial f_{0,i}}{\partial E}}{\omega_{t,i}^2-\omega^2}\right\rangle_W\bar{\varPhi}. \end{equation}

3.1. Derivation of dispersion relation

To evaluate the integrals over the velocity space given in the general linear GAM dispersion relation equation (2.30), we have to choose an equilibrium distribution function for ions. In this work, we consider it to be a bi-Maxwellian, while a regular Maxwellian is used for electrons. The ion distribution function is normalized such that its integral over velocity space equals one ($n_{0,i}=1$). We take $J_{0,i}=1$ (drift kinetic limit):

(3.1)\begin{equation} f_{0,i}=\left(\dfrac{m_i}{2{\rm \pi} T_{\|}}\right)^{{1}/{2}}\left(\dfrac{m_i}{2{\rm \pi} T_\perp}\right)\exp{\left[-\dfrac{m_i}{2}\left(\dfrac{v_{\|}^2}{T_{\|}} +\dfrac{v_\perp^2}{T_\perp}\right)\right]}. \end{equation}

By defining an equivalent temperature, $T_i={T_{\|,i}}/{3}+{2T_{\perp,i}}/{3}$, the bi-Maxwellian can be written in the form

(3.2)\begin{equation} f_{0,i}=\frac{ b^{{3}/{2}}}{{\rm \pi}^{{3}/{2}}v_t^3\chi} \exp{\left[{-}b\left(\frac{v_{\|}^2+v_\perp^2\chi^{{-}1}}{v_t^2}\right)\right]}, \end{equation}

where $b=({2\chi +1})/{3}$, with $\chi ={T_{\perp,i}}/{T_{\|,i}}$, $v_t=\sqrt {{2T_i}/{m}}$ and $E=({m}/{2})(v_{\|}^2+v_\perp ^2)$. We have

(3.3)\begin{equation} \frac{\partial f_{0,i}}{\partial E}={-}\frac{b}{T_i}f_{0,i}. \end{equation}

Equation (2.26) then reduces to

(3.4)\begin{equation} \tilde{\varPhi}_s=\dfrac{\omega\left\langle \dfrac{\bar{\omega}_{d,i}f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W}{1+\dfrac{1}{\tau b}+\omega^2\left\langle \dfrac{f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W}\bar{\varPhi}, \end{equation}

where $\tau ={T_e}/{T_i}$. Substituting this result in the general dispersion relation, we have

(3.5)\begin{equation} \dfrac{2m_ic^2k_r^2}{B^2}+\dfrac{e^2b}{T_i}\left[\left\langle \dfrac{\bar{\omega}_{d,i}^2f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W-\dfrac{\omega^2\left\langle \dfrac{\bar{\omega}_{d,i}f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W^2}{1+\dfrac{1}{\tau b}+\omega^2\left\langle \dfrac{f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W}\right]=0. \end{equation}

These three velocity integrals once evaluated read

(3.6)\begin{gather} \left\langle \dfrac{\bar{\omega}_{d,i}f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W =\dfrac{cm_ik_rv_t^2}{eB_0R_0b\omega_0^2y}\left[y+\left(\dfrac{\chi}{2}+y^2\right)Z(y)\right], \end{gather}

(3.7)\begin{gather} \left\langle \dfrac{\bar{\omega}_{d,i}^2f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W =\left(\dfrac{cm_ik_rv_t^2}{eB_0R_0b}\right)^2\dfrac{1}{\omega_0^2y}\left[\dfrac{y}{2}+y^3+y\chi+\left(\dfrac{\chi^2}{2}+\chi y^2+y^4\right)Z(y)\right], \end{gather}

(3.8)\begin{gather} \left\langle \dfrac{f_{0,i}}{\omega_t^2-\omega^2}\right\rangle_W =\dfrac{1}{\omega_0^2}\dfrac{Z(y)}{y}, \end{gather}

where $Z(y)$ is the plasma dispersion function and

(3.9a,b)\begin{equation} y=\dfrac{\omega}{\omega_0}, \quad \omega_0=\dfrac{v_t}{qR_0\sqrt{b}}. \end{equation}

Then (3.5) becomes

(3.10)\begin{equation} y+q^2\left[F(y)-\dfrac{N^2(y)}{D(y)}\right]=0 , \end{equation}

with

(3.11)\begin{gather} F(y)=\dfrac{y}{2}+y^3+y\chi+\left(\dfrac{\chi^2}{2}+y^2\chi+y^4\right)Z(y), \end{gather}

(3.12)\begin{gather} N(y)=y+\left(\dfrac{\chi}{2}+y^2\right)Z(y), \end{gather}

(3.13)\begin{gather} D(y)=\dfrac{1}{y}\left(1+\dfrac{1}{\tau b}\right)+Z(y). \end{gather}

3.2. Comparison with fluid limit and with previous results

In the limit $\chi =1$, we recover the GAM dispersion relations in Zonca & Chen (Reference Zonca and Chen2008), Zonca et al. (Reference Zonca, Chen and Santoro1996) and Girardo et al. (Reference Girardo, Zarzoso, Dumont, Garbet, Sarazin and Sharapov2014). In the fluid limit, we can show that the dispersion relation equation (3.10) reduces to

(3.14)\begin{equation} \omega^2=\left(\frac{3}{4}+\frac{\chi}{2}+\frac{\chi^2}{2}+b\tau\left(\frac{\chi^2}{4} +\frac{\chi}{2}+\frac{1}{4}\right)\right)\frac{v^2_t}{bR^2_0}. \end{equation}

If we consider the isotropic limit of this fluid dispersion relation, we obtain $\omega ^2_{GAM}=({7}/{4}+\tau )({v^2_t}/{R^2_0})$, which is the same GAM dispersion relation obtained using magnetohydrodynamics with a double adiabatic closure (Smolyakov et al. Reference Smolyakov, Garbet, Falchetto and Ottaviani2008). The GAMs are special types of sound waves and, as sound waves, their frequency strongly depends on the equilibrium pressure. Considering the magnetohydrodynamic description of GAMs with a double adiabatic closure, $\chi$ in our work is equivalent to the ratio ${p_{\perp,0}}/{p_{\|,0}}$, with $p_{\perp,0}$ the equilibrium perpendicular pressure and $p_{\|,0}$ the equilibrium parallel pressure. So increasing $\chi$ for a fixed $p_{\|,0}$ is equivalent to changing the equilibrium perpendicular pressure, which directly modifies the GAM frequency due to its dependence on the equilibrium pressure. This explains the global growing dependence of both the growth rate and frequency on $\chi$, which is evident in the figures presented below. The GAM frequency can be written in terms of the equilibrium parallel and perpendicular pressures as follows (neglecting $\tau$):

(3.15)\begin{equation} \omega^2_{GAM}=\left(\frac{3}{2}+\frac{p_{{\perp},0}}{p_{\|,0}} +\frac{p^2_{{\perp},0}}{p^2_{\|,0}}\right)\frac{p_{\|,0}}{\rho_0 R^2_0}, \end{equation}

where $\rho _0$ is the equilibrium mass density. The dispersion relation equation (3.10) can be written in the following form by substituting (3.11), (3.12) and (3.13) in (3.10):

(3.16)\begin{gather} \frac{1}{b\tau}\left[\frac{1}{q^2}+\frac{1}{2}+\chi+y^2+\left(y^3+y\chi+\frac{\chi^2}{2y}\right)Z(y)\right]\nonumber\\ +\left[\frac{1}{q^2}+\frac{1}{2}+\chi+\frac{yZ(y)}{q^2}+\left(\frac{y}{2}+y\chi+\frac{\chi^2}{2y}\right)+\frac{\chi^2}{4}{Z(y)}^2\right]=0. \end{gather}

3.3. Asymptotic behaviour: case $\tau$ $\rightarrow$ $0$

In this limit, the first term in the square brackets is large compared with the second term. Neglecting this second term, we recover the dispersion relation in Ren & Cao (Reference Ren and Cao2014):

(3.17)\begin{equation} \frac{1}{q^2}+\frac{1}{2}+\chi+y^2+\left(y^3+y\chi+\frac{\chi^2}{2y}\right)Z(y)=0. \end{equation}

Figure 1 shows the frequency and growth rate obtained from the complete dispersion relation (3.10) (blue curve) and the same quantities obtained in the $\tau \rightarrow 0$ limit (red curve).

Figure 1.

(a) Frequency. (b) Growth rate.

3.4. Asymptotic behaviour: case $\tau \rightarrow \infty$

In this limit, the terms in the second square brackets in (3.16) are larger in comparison with those in the first. So the dispersion relation reduces to

(3.18)\begin{equation} \frac{1}{q^2}+\frac{1}{2}+\chi+\frac{yZ(y)}{q^2}+\left(\frac{y}{2}+y\chi+\frac{\chi^2}{2y}\right)+\frac{\chi^2}{4}{Z(y)}^2=0. \end{equation}

Figure 2 shows the frequency and growth rate obtained from the complete dispersion relation (3.10) (blue curve) and those obtained in the $\tau \rightarrow \infty$ limit (red dashed curve).

Figure 2.

(a) Frequency. (b) Growth rate.

In the following sections, we solve the complete dispersion relation, (3.10), with values of $\tau$ that are more realistic for tokamak plasmas.

3.5. Effect of ion temperature anisotropy on GAM frequency and growth rate

The GAM frequency and growth rate increase with $\chi$ (figure 3). However, the growth rate saturates for values of the parameters for which the instability becomes marginally stable. This saturation occurs at lower values of $\chi$ for higher values of $q$. It should be noted that the growth rate at higher values of $q$ is overestimated, since damping effects due to finite orbit width are not considered in our model. These effects tend to be more important when $q$ increases (Biancalani et al. Reference Biancalani, Bottino, Lauber and Zarzoso2014). Similar results were obtained in Ren & Cao (Reference Ren and Cao2014). The regime of validity of our model is the plasma core. This is the region of excitation of energetic-particle-induced GAMs (EGAMs). Incidentally, the zonal flow residual level (Cho & Hahm Reference Cho and Hahm2021) can be enhanced by the anisotropic energetic particle distribution dominated by barely passing/barely trapped particles or by deeply trapped particles (Lu et al. Reference Lu, Wang, Lauber, Fable, Bottino, Hornsby, Hayward-Schneider, Zonca and Angioni2019).

Figure 3.

(a) Frequency. (b) Growth rate.

3.6. Effects of electron to ion temperature ratio

In this section, we study the effect of a finite $\tau$ (this parameter was neglected in Ren & Cao Reference Ren and Cao2014). Figure 4(a,b) shows the effects of $\chi$ on GAM frequency and growth rate for two different values of $\tau$. We observe, as expected, that these quantities are increasing functions of $\chi$. However, there is a significant increase in frequency and growth rate with $\tau$. Figure 4(c,d) shows the variation of GAM frequency and growth rate with $\tau$ for two values of $\chi$. We recover the isotropic $\tau$ dependence of GAMs (i.e. increasing frequency with increasing $\tau$). This effect is more pronounced for higher values of $\chi$ in the anisotropic case. Figure 5 shows the effect of $\tau$ on GAM frequency and growth rate for two different values of the safety factor $q$ and for a fixed $\chi$. We observe stronger damping for smaller values of $q$. This confirms the fact that even in the anisotropic case, GAMs are more stable in the core than at the edge of the tokamak plasma, where the safety factor $q$ has large values. We can conclude from figures 4 and 5 that GAM frequency and growth rate significantly increase with $\tau$ and neglecting this parameter can lead to an underestimation of the frequency and to an overestimation of the damping rate.

Figure 4.

Effects of $\tau$: (a,c) frequency; (b,d) growth rate.

Figure 5.

Effects of $\tau$: (a) frequency; (b) growth rate.