2.1. A way to derive eigenmode equations in warm plasmas; reconciling standing and travelling waves

Eigenmodes in bounded plasmas are standing waves, but these waves represent a superposition of travelling waves. This statement is employed below for both warm and cold plasmas.

As known, eigenmodes can be described by the following equation (see e.g. Stix (Reference Stix1992)):

(2.1)\begin{equation} \boldsymbol{\nabla}\times\boldsymbol{\nabla} \times \delta \boldsymbol E-\frac{\omega^2}{c^2}{ \stackrel{\leftrightarrow}{\boldsymbol\epsilon} } \delta \boldsymbol E=0, \end{equation}

where $\delta \boldsymbol E\propto \exp (-{\rm i}\omega t)$ is a perturbed electric field, $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ is a dielectric tensor. This equation immediately follows from the Maxwell equations $\boldsymbol {\nabla }\times \delta \boldsymbol E=({{\rm i}\omega }/{c})\delta \boldsymbol B$ and $\boldsymbol {\nabla }\times \delta \boldsymbol B=-({{\rm i}\omega }/{c})\stackrel {\leftrightarrow }{\boldsymbol\epsilon }\delta \boldsymbol E$, where $\delta \boldsymbol B$ is the perturbed magnetic field. Equation (2.1) should be supplemented by boundary conditions. Then it will determine waves standing in the radial direction (in cylindrical geometry). When $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ does not depend on the radial structure of perturbation, as in cold plasmas, (2.1) reduces to a second-order differential equation.

However, in general this is not the case: $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ is determined by the disturbed distribution function of the particles, $\delta f$, which, in turn, is determined by $\delta \boldsymbol {E}$,

(2.2)\begin{equation} \delta f =-\frac{e}{M}\int_{-\infty}^t\,{\rm d}t^\prime \left [\delta \boldsymbol{E}- \frac{{\rm i}}{\omega} \boldsymbol{v}\times\boldsymbol{\nabla} \times \delta \boldsymbol{E} \right ]\frac{\partial F }{\partial \boldsymbol{v}}, \end{equation}

where $F$ is an equilibrium distribution function, $e$ and $M$ are the particle charge and mass, respectively, $\boldsymbol r =\boldsymbol r (t^\prime )$ and $\boldsymbol {v} =\boldsymbol {v} (t^\prime )$ (the integral is taken along particle orbits). The known technique for the calculation of integral in (2.2) with Maxwell distribution function ($F_M$) relies on perturbations in the form of a plane travelling wave,

(2.3)\begin{equation} \delta\boldsymbol E\propto \exp({\rm i}\boldsymbol{k}\boldsymbol{r} -{\rm i}\omega t). \end{equation}

This means that $\delta \boldsymbol {E}\propto {\rm e}^{{\rm i}k_rr}$, which is not consistent with $\delta \boldsymbol {E}(r)$ determined by (2.1) ($r$ is the radial coordinate). However, the inconsistency diminishes when $k_r\ll k_\vartheta$ ($\vartheta$ is the poloidal angle): the dependence of $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ on $k_r$ becomes negligible, which implies that $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ is valid for various $k_r$, not only for the radial wavenumber of (2.3). To see this, we note that the propagation of the wave (2.3) across the magnetic field manifests itself in $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ in the form of Bessel functions with the argument that contains Larmor radius in the product $k_\perp \rho _T$ ($\rho _T$ is the Larmor radius of thermal particles and $k_\perp$ is transverse wavenumber). Due to this, $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ is valid for any $k_r$ that satisfies the condition $k_r\ll k_\vartheta$. On the other hand, when $k_\perp \rho _T \rightarrow 0$ due to vanishing Larmor radius, the condition $k_r\ll k_\vartheta$ is not necessary. Thus, (2.1) with $\stackrel {\leftrightarrow }{\boldsymbol\epsilon }$ calculated for a plane wave is self-consistent when $\rho _T$ is neglected (in cold plasmas); it is self-consistent for $\rho _T\neq 0$ provided that $k_\perp ^2\approx k_\vartheta ^2$. Otherwise, there is a conflict between (2.1) and (2.3). Nevertheless, there is a way to use the technique developed for plane travelling waves in every case: for this purpose, one has to expand $\delta \boldsymbol {E}$ in a Fourier series with terms in the form of plane waves. This approach will be used in this work.

Below we consider an example which demonstrates that the mentioned Fourier expansion can be useful even in the analysis of eigenmodes in a cold plasma.

The studies of edge-localized ICE associated with FMM often relied on a simple model employing plasma inhomogeneity to provide a potential well in a Schrödinger-like equation, see Coppi (Reference Coppi1993) and the overview Gorelenkov (Reference Gorelenkov2016). In this model the eigenmode amplitude is a Gaussian,

(2.4)\begin{equation} \delta{E}(x) =\delta{E}_0\exp\left [-\frac{(x-x_0)^2}{\varDelta ^2}\right ], \end{equation}

where $x=r/a$, $a$ is the plasma radius, and $x_0$ and $\varDelta$ are parameters defined by equations (7)–(9) in Coppi (Reference Coppi1993). It was obtained in the assumption that the radial derivative terms in the eigenmode equation are small, $\partial /\partial r \ll m/r$ ($m$ is the poloidal mode number). This implies that an effective radial wavenumber defined by $ik_r^{ef}\delta E\sim \partial \delta E/\partial r$ is small ($k_r^{ef}\ll k$) and, therefore, the FMM frequency for small longitudinal wavenumber ($k_\|$) is approximately determined by $k_\vartheta =m/r$ (Coppi Reference Coppi1993),

(2.5)\begin{equation} \omega^2 = k_\vartheta^2(x_0)v_A^2(x_0)+\delta\omega^2, \end{equation}

where $\delta \omega ^2\ll \omega ^2$. However, the exact magnitude of $k_{r}^{ef}$ cannot be introduced and, therefore, it is not clear how large $m$ should be for (2.4), (2.5) to hold.

In order to clarify this point, let us write $\delta {E}(x)$ as a superposition of plane waves with the radial wavenumbers $k_{r,\nu }=2{\rm \pi} \nu /a$, where $\nu =0, 1, 2,\ldots,\infty$:

(2.6)\begin{equation} \delta{E}(x) =\frac{A_0}{2} +\sum_{\nu=1}^{\infty} \left [A_\nu\cos (2{\rm \pi} \nu x) +B_\nu\sin (2{\rm \pi} \nu x) \right ], \end{equation}

where

(2.7a,b)\begin{equation} A_\nu =2\int_0^1 \,{{\rm d} x} \delta E(x)\cos (2{\rm \pi} \nu x),\quad B_\nu =2\int_0^1 \,{{\rm d} x} \delta E(x)\sin (2{\rm \pi} \nu x). \end{equation}

The infinite series in (2.6) is consistent with the condition $k_r\ll k_\vartheta$ only when the terms with large $\nu$ weakly contribute and can be neglected. To consider this case, we replace $\sum _0^{\infty }$ with $\sum _0^{\nu _{\max }}$, assuming that terms with $\nu >\nu _{\max }$ are negligible. Then $k_{r,\nu } \leq k_{r,\max }=2{\rm \pi} \nu _{\max } /a$, and the condition $k_{r,\max }\ll k_\vartheta$ requires

(2.8)\begin{equation} \nu_{\max}\ll \frac{m}{2{\rm \pi} x}. \end{equation}

Because $\nu _{\max }=\mathcal {N}-1$ ($\mathcal {N}$ is the number of the terms in the sum with $\nu \leq \nu _{\max }$), (2.8) represents the mildest restriction on $m$ when $\mathcal {N}$ is smallest. Thus, we need to know the minimum number $\mathcal {N}_{\min }$ of the terms in (2.6) which provides a good approximation for (2.4). After $\mathcal {N}_{\min }$ is found, by using (2.8) we can obtain the poloidal mode numbers for which $k_r\ll k_\vartheta$:

(2.9)\begin{equation} m\gg m_{{\rm crit}} \equiv 2{\rm \pi} x \left(\mathcal{N}_{\min} -1\right). \end{equation}

As in Coppi (Reference Coppi1993), we consider the JET preliminary tritium experiment where the ICE was observed (Cottrell et al. Reference Cottrell, Bhatnagar, Costa, Dendy, Jacquinot, McClements, McCune, Nave, Smeulders and Start1993), and take $x_0=r_0/a=0.816$, $\varDelta =0.738/\sqrt {m}$, which correspond to ion density profile $n_i=n_i(0)(1-x^2)^{1/2}$. We found that the mode radial structure (2.4) can be approximated satisfactory by a superposition of plane waves with $k_{r,\nu }\ll k_\vartheta$ only for $m\gg m_{{\rm crit}} =5(\mathcal {N}_{\min } -1)\gtrsim 15$, with $\mathcal {N}_{\min }\sim 4\unicode{x2013}5$, see figure 1. Taking $\omega \approx m v_A/r\approx l\omega _{Bi}$ (with $r=0.85$ m, $v_A =10^7\,{\rm m}\,{\rm s}^{-1}$, the ion gyrofrequency $\omega _{Bi}=1.35\times 10^8$ s$^{-1}$ in deuterium plasma, and $l$ is an integer) and $m\gg m_{{\rm crit}}$ we obtain $l\gg 1$. This means that the model of Coppi (Reference Coppi1993) is not suitable for the description of cyclotron modes with $l\sim 1$ for the parameters used.

Figure 1. Approximation of the amplitude of the mode given by (2.4) with $m=50$ (black line with markers) by several Fourier coefficients (colour lines). We observe that five harmonics ($\nu _{\max }=4$) provide a good approximation.

Note that although all works employing the described model (and its modifications by including effects of toroidicity, etc.) assumed $m\gg 1$, no attempt to find $m_{{\rm crit}}$ was done yet. On the other hand, it was found recently that FMMs with frequencies $\omega \geq \omega _{Bi}$, including $\omega \approx l\omega _{Bi}$, can occupy a considerable part of the plasma cross-section and can have maximum amplitudes in the plasma core (Burdo & Kolesnichenko Reference Burdo and Kolesnichenko2020). The restriction $m\gg m_{{\rm crit}}$ is absent for these modes. However, the analysis of the mentioned work is not applicable to warm plasmas.

2.2. Toroidicity parameters and choice of approximation

When the magnetic field is homogeneous, MHD modes in a cold plasma – Alfvén modes with frequency $\omega =k_\| v_A$ and FMM with $\omega =k v_A$ – can be obtained in a kinetic theory provided that

(2.10)\begin{gather} x_l\equiv\frac{\omega -l\omega_{B,j}}{k_\|v_{T,j}}\gg 1, \end{gather}

(2.11)\begin{gather} k_\perp\rho_{T,j} \ll 1, \end{gather}

where $v_T =\sqrt {2T/M}$ is thermal velocity, $\rho _T =v_T/\omega _B$, the subscripts $\|$ and $\perp$ label magnitudes along and across the magnetic field, respectively, $j =e,i$ labels electrons and ions. Wave damping is exponentially small or vanishing due to the first inequality. The Larmor radius can be eliminated from $\epsilon _{ij}$ when the second inequality holds. Finite (non-vanishing) $k_\perp \rho _T$ leads to the existence of modes that are absent in MHD and in the two-fluid model. Because the Larmor radius $\rho _\perp =v_\perp /\omega _B$ ($v$ is the particle velocity) is responsible for the deflection of individual particles from the magnetic field lines, kinetic theory contains only one parameter associated with transverse motion,

(2.12)\begin{equation} \xi=k_\perp\rho_\perp. \end{equation}

The inhomogeneity of the magnetic field increases the number of characteristic parameters, as described below.

In tokamaks, the magnetic field can be approximated by

(2.13)\begin{equation} B=\bar{B}(1-\epsilon \cos \vartheta), \end{equation}

where $\bar B$ is the magnetic field on the magnetic axis, $\vartheta$ is the poloidal angle, $\epsilon =r/R$, $R$ is the major radius of the torus. In this field, the particle radial excursion is $\varDelta r\sim v_D\tau _b$, where $v_D\sim v\rho /R$ is the toroidal drift velocity, $\rho =v/\omega _B$, $\tau _b\sim \omega _b^{-1}$, $\omega _b\sim v_\|/(qR)$ is transit/bounce frequency, $q$ is the tokamak safety factor. This leads to the following characteristic drift parameter for passing and trapped particles:

(2.14)\begin{equation} \xi_D= k_\perp v_D\tau_b. \end{equation}

In addition, bounce oscillations of trapped particles along the magnetic field result in

(2.15)\begin{equation} \xi_{b}= k_\|\mathcal{A}, \end{equation}

where $\mathcal {A}\sim v_\| \tau _b \sim 2\kappa qR$ is the oscillation amplitude, $2\kappa$ represents turning points $\vartheta _t$ of trapped particles, $\kappa =[\mathcal {E} -\mu \bar {B}(1-\epsilon )]/(2\mu \bar {B}\epsilon )$ is the particle trapping parameter ($\kappa <1$ for trapped particles), $\mathcal {E}$ and $\mu$ are the particle energy and magnetic moment, respectively. One more parameter arises because of the dependence of $\omega _B$ on $\vartheta$. For passing particles moving along the field lines this parameter is

(2.16)\begin{equation} \zeta_l = \epsilon l\bar\omega_B\tau_b=lq\frac{r}{\rho_\|}, \end{equation}

where $\rho _\|=v_\|/\bar \omega _B$.

All these parameters are a consequence of the periodic motion of the particles. They appear in arguments of Bessel functions as a result of the expansion of the electromagnetic field perturbations proportional to $\exp (ix\sin \phi )$ in the Fourier series

(2.17)\begin{equation} {\rm e}^{{\pm} {\rm i}x\sin\phi} =\sum_{p=-\infty}^{\infty} {\rm J}_p(x) {\rm e}^{{\pm} {\rm i}p\phi}, \end{equation}

where ${\rm J}_p(x)$ is the Bessel function of the first kind, and $p$ is an integer. For Larmor rotation $\phi =\omega _B t$, $x=\xi$; for bounce/transit motion $\phi =\omega _{b} t$, $x$ is represented by $\xi _D$, $\xi _{b}$ and $\zeta _l$. The left-hand side of (2.17) represents a part of the phase $\varPsi$ of the perturbation $\delta X\propto \exp (-{\rm i}\varPsi )$, with $\varPsi (t) =\int ^t\,{\rm d}t^\prime (\omega -m\dot {\vartheta } +n\dot {\varphi })$, where $\varphi$ is the toroidal angle, a dot over letters denotes a time derivative, $n$ is the toroidal mode number. The right-hand side of (2.17) produces terms proportional to $l\omega _B$ and $s\omega _b$ ($l$ and $s$ are integers) in the phase, which leads to an infinite number of wave–particle resonances.

Although (2.17) is the same for both types of the particle motion, the roles of Larmor rotation and transit/bounce motion are different: transit/bounce motion manifests itself only for long time scale, $\varDelta t\gg \omega _B^{-1}$. Due to this, the approximation of the straight magnetic field is suitable for the study of fast instabilities, i.e. instabilities with the growth rate $\gamma > \omega _b$. Another difference is that many terms of (2.17) are important for transit/bounce motion (for high frequency modes, but not for low frequency modes, like toroidicity-induced Alfvén eigenmodes); in contrast to this, one can consider a separate cyclotron harmonic, at least for $\omega \approx l\omega _{Bi}$ with low $l$ (the resonances of different harmonics overlap at sufficiently large $l$, which can explain the continuum ICE spectrum for $l>7$ (Fülöp et al. Reference Fülöp, Kolesnichenko, Lisak and Anderson1997) observed in JET (Cottrell et al. Reference Cottrell, Bhatnagar, Costa, Dendy, Jacquinot, McClements, McCune, Nave, Smeulders and Start1993)).

Because $\zeta _l$ and $\xi _D$ are relevant to the same time scale, it is of interest to compare them. Taking $\xi _D \sim k_\perp q \rho v/v_\|$ we obtain

(2.18)\begin{equation} \frac{\zeta_l}{\xi_D} \sim \frac{l}{k_\perp \rho} \frac{r}{\rho}. \end{equation}

We conclude that typically $\zeta _l\gg \xi _D$.

The presence of two groups of particles – passing and trapped – considerably complicates eigenmode equations. First, the large time scale $\sim \tau _b$ should lead to two infinite series of Bessel functions. Second, orbits with $\kappa \sim 1$ (marginally trapped and marginally passing orbits) are known to be described by special functions (by elliptic integrals). Third, the presence of a border between trapped and passing particles in the velocity space further complicates equations. These facts make analytical calculation of the integrals in (2.2) impossible, unless simplifying assumptions are made. Note that this problem can be minimized for fast ions: one simplification arises when the fast ion population consists of either trapped or passing particles (not both together). Another simplification occurs when, as it is often the case, the perturbative approach is sufficient. Then one needs to know only the anti-Hermitian part of the fast-ion contribution to the dielectric tensor.

In this work, we restrict ourselves to consideration of passing particles only, neglecting the presence of trapped particles. This in justified when toroidicity is very small (the fraction of trapped particles at a certain flux surface equals $0.9 \sqrt {r/R}$ when orbits are standard and velocity distribution is isotropic). In addition, this can also make sense for any toroidicity provided that $k_\|qR\ll 1$: then $\xi _b \ll 1$ and, therefore, the term with $p=0$ in (2.17) mainly contributes for $x=\xi _b$, which implies that the bounce oscillation parameter is not important. On the other hand, the role of finite orbit width, $\xi _D$, is similar to that for passing particles. The latter is expected to be of minor importance due to (2.18), see § 3.2.

3.1. Equation for eigenmodes when the magnetic field is homogeneous

The approximation of homogeneous magnetic field can be justified by the fact that ICE is a local phenomenon, as shown by contemporary particle-in-cell simulations across most tokamaks and stellarators (Carbajal et al. Reference Carbajal, Dendy, Chapman and Cook2017). Therefore, many results are obtained in this approximation, and we first derive equations for a homogeneous magnetic field.

We proceed from a Poisson equation

(3.1)\begin{equation} \nabla^2\delta{\varPhi} =-\sum_{j=e,i} 4{\rm \pi} e_j\int \,{\rm d}^3v \delta f_j, \end{equation}

where $\delta f$ is given by (2.2) which for $\delta \boldsymbol {B}=0$ reduces to

(3.2)\begin{equation} \delta f_j=\frac{e_j}{M_j}\int_{-\infty}^t \,{\rm d}\tau\boldsymbol{\nabla}\delta\varPhi(\tau)\boldsymbol{\cdot}\frac{\partial F_j}{\partial\boldsymbol{v}(\tau)}. \end{equation}

We take $\delta \varPhi$ in the form of a monochromatic travelling wave in the poloidal and toroidal directions, whereas its radial dependence will remain unspecified but presented as a superposition of plane waves,

(3.3)\begin{gather} \delta\varPhi =\hat{\varPhi} (r) \exp(-{\rm i}\omega t +{\rm i}m\vartheta -{\rm i}n\varphi), \end{gather}

(3.4)\begin{gather} \hat\varPhi (r)=\sum_{\nu =-\infty}^{\infty} \varPhi_\nu {\rm e}^{{\rm i}k_\nu r}, \end{gather}

where $\varPhi _\nu$ are Fourier coefficients defined by

(3.5)\begin{equation} \varPhi_\nu =\frac{1}{a}\int_0^a \,{\rm d}r^\prime \hat{\varPhi}(r^\prime) {\rm e}^{-{\rm i}k_\nu r^\prime}, \end{equation}

and $k_\nu \equiv k_r(\nu )=2{\rm \pi} \nu /a$ are the radial mode numbers of harmonics constituting the mode, and $\nu$ are integers (in contrast to (2.6), here $-\infty <\nu <\infty$). The perturbation (3.4) in the integrand of (3.2) can be written as (we take into account that coordinates and velocities depend on time, as explained after (2.2))

(3.6)\begin{equation} \delta\varPhi (t,\tau) =\sum_\nu \varPhi_\nu \exp({-{\rm i}\omega t+\boldsymbol{k}_\nu\boldsymbol{x}(t)}) \exp\left({-{\rm i}\int_t^\tau \,{\rm d}t^\prime (\omega -\boldsymbol{k}_\nu\boldsymbol{v}(t^\prime))}\right), \end{equation}

where $\boldsymbol {k}_\nu$ is a wavevector with $k_r=k_\nu$, $\boldsymbol {k}_\nu \boldsymbol {v}= k_\nu v_r +k_b v_b +k_\|v_\|$, the subscript ‘$b$’ labels binormal components of the vectors ($k_\| = k_\vartheta b_\vartheta +k_\varphi b_\varphi =(mq^{-1}-n)/R$, $k_b=k_\vartheta b_\varphi -k_\varphi b_\vartheta \approx m/r$ for $m/n \gg \epsilon ^2 /q$, $\boldsymbol b =\boldsymbol B /B$). Then (3.2) for $F=F(\mathcal {E})$ takes the form

(3.7)\begin{equation} \hat{f} = \sum_{\nu}ie\varPhi_\nu \frac{\partial F} {\partial \mathcal{E}} \left [\int_{-\infty}^t \,{\rm d}\tau A_\nu (t,\tau)\boldsymbol{k}_\nu\boldsymbol{v} (\tau)\right ] {\rm e}^{{\rm i}k_\nu r }, \end{equation}

where

(3.8)\begin{equation} A_\nu (t,\tau) = \exp\left [-{\rm i}\int_t^{\tau} \,{\rm d}t^\prime \varOmega (t^\prime)\right ], \end{equation}

$\varOmega =\omega -\boldsymbol {k}_\nu \boldsymbol {v}$, the subscript $j$, which labels species, is omitted. Because $\boldsymbol {k}_\nu \boldsymbol {v} =\omega -\varOmega$ and $\varOmega (\tau )A_\nu (t,\tau )={\rm i}\partial A_\nu (t,\tau )/\partial \tau$, we can integrate the term proportional to $\varOmega$ in (3.7). Then, due to $A_\nu (t,t) =1$,

(3.9)\begin{equation} \hat{f} = \sum_{\nu}e\varPhi_\nu\frac{\partial F} {\partial \mathcal{E}} \left [{\rm i}\omega\int_{-\infty}^t \,{\rm d}\tau A_\nu (t,\tau) +1\right ]{\rm e}^{{\rm i}k_\nu r }. \end{equation}

In order to calculate the integral $\int _t^\tau \boldsymbol {k}_\nu \boldsymbol {v}\,{\rm d}t^\prime$ in (3.8) we write $\boldsymbol {k}_{\perp,\nu }\boldsymbol {v}_\perp =k_{\perp,\nu }v_\perp \cos (\alpha -\psi _\nu )$, with $k_{\perp,\nu }^2=k_\nu ^2+k_b^2$, $\psi _\nu$ the angle defined by $k_\nu =k_{\perp,\nu }\cos \psi _\nu$, and $\alpha$ the gyroangle, $\dot \alpha =-\omega _B$ (see e.g. Belikov & Kolesnichenko Reference Belikov and Kolesnichenko1982; Mikhailovskii Reference Mikhailovskii1986). This yields

(3.10)\begin{equation} \int_t^\tau \,{\rm d}t^\prime k_{{\perp},\nu} v_\perp\cos (\alpha -\psi)= \xi_\nu(\sin \alpha_t -\sin \alpha_\tau), \end{equation}

with $\xi _\nu =k_{\perp,\nu } v_\perp /\omega _B$, $\sin \alpha _t \equiv \sin (\alpha (t)-\psi )$ and $\sin \alpha _\tau \equiv \sin (\alpha (\tau )-\psi )$. Here $\sin \alpha _\tau$ can be transformed to $\alpha _\tau$ by Fourier expansion (2.17), with

(3.11)\begin{equation} \alpha_\tau-\alpha_t=\int_t^\tau \,{\rm d}t^\prime \dot\alpha =\int_\tau^t \,{\rm d}t^\prime\omega_B. \end{equation}

As a result, we obtain

(3.12)\begin{equation} A_\nu (t,\tau) = \sum_l {\rm J}_l(\xi_\nu)\exp({{\rm i}\xi_\nu\sin\alpha_t-{\rm i}l\alpha_t}) \exp\left({{\rm i}\int_\tau^t \,{\rm d}t^\prime\varOmega_l}\right), \end{equation}

where $\varOmega _l =\omega -l\omega _B -k_\|v_\|$. Because $\varOmega _l= {\rm const.}$ in (3.12), $A(t,\tau ) =A(t-\tau )$, which leads to

(3.13)\begin{equation} \int_{-\infty}^t \,{\rm d}\tau A_\nu (t,\tau) = {\rm i}\sum_l\mathcal{Z} (\varOmega_l) {\rm J}_l(\xi_\nu) \exp({{\rm i}\xi_\nu\sin\alpha_t-{\rm i}l\alpha_t}), \end{equation}

where

(3.14)\begin{equation} \mathcal{Z} (\varOmega_l)=-{\rm i}\int_0^\infty {\rm e}^{{\rm i}\varOmega_lt^\prime} \,{\rm d}t^\prime =\frac{P}{\varOmega_l}-{\rm i}{\rm \pi} \delta (\varOmega_l), \end{equation}

and $P/\varOmega _l$ means Cauchy principal value. Combining (3.13), (3.9) and using relation

(3.15)\begin{equation} \oint \frac{{\rm d}\phi}{2{\rm \pi}} \exp({\pm {\rm i}x\sin \phi \mp {\rm i}n\phi}) ={\rm J}_n(x), \end{equation}

we obtain after alpha averaging

(3.16)\begin{equation} \oint \frac{{\rm d}\alpha}{2{\rm \pi}}\hat f=-\sum_{l,\nu} e{\varPhi}_\nu\frac{{\rm d}F}{{\rm d}\mathcal{E}} \left [\mathcal{Z} (\varOmega_l){\rm J}_l^2(\xi)\omega -1 \right]. \end{equation}

For Maxwell distribution $F_M=n(r)F_\|(v_\|)F_\perp (v_\perp )$ ($n(r)$ is the particle density) the known integrals are

(3.17)\begin{equation} 2{\rm \pi} \int_0^\infty \,{\rm d}v_\perp v_\perp F_\perp {\rm J}_l^2(\xi)=\hat{{\rm I}}_l(z) \end{equation}

and

(3.18)\begin{equation} \int_{-\infty}^{\infty}\,{\rm d}v_\|\mathcal{Z}(\varOmega_l)F_\|= \frac{1}{\omega-l\omega_B}Z_l(x_l). \end{equation}

Here

(3.19)\begin{equation} Z_l(x_l)=\int_{-\infty}^{\infty}\,{\rm d}v_\|\frac{F_\|}{1-v_\|/(v_Tx_l)}, \end{equation}

$\hat {{\rm I}}_l(z) \equiv {\rm I}_l(z){\rm e}^{-z}$, ${\rm I}_l(z)$ is the modified Bessel function of the first kind, $z=0.5k_\perp ^2 \rho _T^2$, $x_l$ is given by (2.10); $Z_l(x_l)\approx 1$ for $x_l\gg 1$ and $Z_l(x_l)\approx 2x_l^2-{\rm i}\sqrt {{\rm \pi} }x_l$ for $x_l\ll 1$ (Shafranov Reference Shafranov1967). Note that $Z_l(x) =-xZ_{dis}(x)$, where $Z_{dis}(x)$ is the plasma dispersion function of Richardson (Reference Richardson2019).

Now we are ready to finalize the eigenmode equation. We write (3.1) as

(3.20)\begin{equation} \frac{1}{r}\frac{{\rm d}}{{\rm d}r}r\frac{{\rm d} \hat{\varPhi}}{{\rm d} r} - \left (\frac{m^2}{r^2} + \frac{n^2}{R^2}\right )\hat{\varPhi} =-\sum_{j=e,i} 4{\rm \pi} e_j\int\, {\rm d}\boldsymbol{v}_\perp \,{\rm d}v_\| \hat{f}_j, \end{equation}

where ${\rm d}\boldsymbol {v}_\perp ={\rm d}v_\perp v_\perp \,{\rm d}\alpha$. Using (3.16)–(3.20) we obtain

(3.21)\begin{align} & \frac{1}{r}\frac{{\rm d}}{{\rm d}r}r\frac{{\rm d} \hat{\varPhi}}{{\rm d} r} - \left (\frac{m^2}{r^2} + \frac{n^2}{R^2}\right )\hat{\varPhi} \nonumber\\ & \qquad +\sum_{j=e,i}\sum_{\nu=-\infty}^{\infty}\frac{1}{d_j^2}\left [\sum_{l=-\infty}^{\infty}\hat{{\rm I}}_l(z_{\nu,j})Z_l(x_l)\frac{\omega}{\omega -l\omega_B} -1\right ]\varPhi_\nu {\rm e}^{{\rm i}k_\nu r}=0, \end{align}

where $\varPhi _\nu$ is given by (3.5). Normally a finite number of $\nu$-harmonics is important (see the example in § 2.1). Therefore, when $z_{\nu =1,e}\ll 1$, $z_{\nu \neq 1,e}$ can be small, too. Then the electron contribution to the mode equation reduces to

(3.22)\begin{equation} \frac{1}{d_e^2}\left [Z_0(x_0) -1 \right ]\hat{\varPhi}. \end{equation}

Equation (3.21) can be transformed to an equation for a plane travelling wave by keeping only one term in the sum over $\nu$ and replacing the first term with $-k_r^2$. In this particular case (3.21) coincides with the equation which can be obtained from that of Harris (Reference Harris1961) after applying relations (3.17), (3.18).

3.2. Equation with a non-uniform particle gyration

One can expect that the main effect of toroidicity will be associated with a non-uniform gyration of passing particles: as shown in § 2.2, $\zeta _l \gg \xi _D$, unless $k_\perp \rho$ is unrealistically large. Below we consider this case.

We employ (3.12) with $\xi =\bar \xi$ (bar over letters means that $B=\bar {B}$), which can be justified as follows.

As in a homogeneous magnetic field, the gyration can be described by (Mikhailovskii Reference Mikhailovskii1986; Belikov & Kolesnichenko Reference Belikov and Kolesnichenko1982)

(3.23)\begin{equation} \tilde{r}(\tau) - \tilde{r}(t)\approx \frac{v_\perp}{\omega_B}\sin\alpha (t)-\frac{v_\perp}{\omega_B}\sin\alpha (\tau). \end{equation}

This agrees with the right-hand side of (3.10) and, therefore, we can use (3.10) even for inhomogeneous magnetic field. However, now we have to take into account that $\alpha (t)$ is not a linear function ($\dot \alpha \neq$ const.). The variation of $\dot \alpha$ affects the phase in (3.10) for a particle transit time considerably. This follows from (3.11) and the relations

(3.24a,b)\begin{equation} \varOmega_l=\bar{\varOmega}_l-l\bar{\omega}_B\epsilon \cos\vartheta,\quad \bar{\varOmega}_l=\omega-l\bar{\omega}_B -k_\|v_\|. \end{equation}

There is one more effect of the inhomogeneity of the magnetic field: the transverse Larmor radius $\rho _\perp \equiv v_\perp /\omega _B$ varies when $B=B(\vartheta )$. Using conservation of the particle magnetic moment we find that $\rho _\perp =\bar \rho _\perp (1+0.5\epsilon \cos \vartheta )$. Due to this, (3.10) leads to (3.12) with $\xi =\bar \xi$ provided $0.5\epsilon \bar \xi$ is sufficiently small, so that $\exp [{\rm i}\epsilon \bar \xi /2] \approx 1$. Below we assume that this is the case and omit the bar over $\xi$.

To calculate the integral in (3.12) we note that

(3.25)\begin{equation} \exp({{\rm i}\int_{\tau}^t \,{\rm d}t^\prime \varOmega_l (t^\prime)})= \exp({{\rm i}\bar{\varOmega}_l(t-\tau) -{\rm i}\zeta_l\sin \vartheta (t)}) \exp({{\rm i}\zeta_l\sin \vartheta (\tau)}), \end{equation}

where $\zeta _l =l\bar {\omega }_B qr/v_\|$ and $\sin \vartheta (\tau )$ can be transformed to $\tau$ by means of the Bessel series expansion (2.17). Then due to the relation

(3.26)\begin{equation} \vartheta (t) -\frac{v_\|}{qR}t = \vartheta (\tau) -\frac{v_\|}{qR}\tau, \end{equation}

which results from the equation of motion for the particle guiding centre, we obtain

(3.27)\begin{equation} \int_{-\infty}^t \,{\rm d}\tau \exp\left({{\rm i}\int_{\tau}^t \,{\rm d}t^\prime\varOmega_l(t^\prime)}\right)= {\rm i}\sum_s\frac{{\rm J}_s(\zeta_l)}{\bar{\varOmega}_l-sv_\|/(qR)} \exp\left({-{\rm i}\zeta_l\sin\vartheta (t)+{\rm i}s\vartheta (t)}\right). \end{equation}

The flux surface averaging and the gyrophase averaging result in

(3.28)\begin{equation} \oint\frac{{\rm d}\vartheta}{2{\rm \pi}}\oint\frac{{\rm d}\alpha}{2{\rm \pi}}\int_{-\infty}^t A_\nu(t,\tau) ={\rm i}\sum_{l,s}\frac{{\rm J}_l^2(\xi){\rm J}_s^2(\zeta_l)}{\bar{\varOmega}_{l,s}}, \end{equation}

where

(3.29)\begin{equation} \bar{\varOmega}_{l,s}=\omega-l\bar{\omega}_B-k_{\|,s}v_\|, \end{equation}

with $k_{\|,s}=k_\|+s/(qR)$.

For Maxwell distribution, the velocity integral of (3.28) can be written as

(3.30)\begin{equation} \int \,{\rm d}\boldsymbol{v}_\perp F_\perp \int \,{\rm d}v_\parallel F_\| \oint \frac{{\rm d}\vartheta}{2{\rm \pi}} \oint \frac{{\rm d}\alpha}{2{\rm \pi}} \int_{-\infty}^t \,{\rm d}\tau A_\nu (t,\tau)={\rm i}\sum_l {\rm Y}_l\hat{{\rm I}}_l(z_\nu), \end{equation}

where

(3.31)\begin{equation} {\rm Y}_l=\sum_{s=-\infty}^\infty\int \,{\rm d}v_\| F_\| \frac{{\rm J}_s^2(\zeta_l)}{\bar{\varOmega}_{l,s}}. \end{equation}

In the limit case of a plasma with zero temperature $F_\|=\delta (v_\|)$ ($\delta (x)$ is a Dirac delta function) and, therefore, $\bar \varOmega _{l,s}$ is not involved in the summation. Due to this, taking into account that

(3.32)\begin{equation} \sum_{s=-\infty}^{\infty} {\rm J}_s^2 =1, \end{equation}

we obtain

(3.33)\begin{equation} {\rm Y}_l(T=0) = \frac{1}{\omega -l\bar\omega_B}. \end{equation}

The same relation takes place in homogeneous magnetic field: $\int \,{\rm d}v_\| F_\|(T=0)/(\omega -l\bar \omega _B-k_\|v_\|)=1/(\omega -l\bar \omega _B)$. This is explained by the fact that there is no particle motion along the field lines in the cold plasma.

It is convenient to introduce $\mathcal {Y}_l$ by (compare with (3.18))

(3.34)\begin{equation} {\rm Y}_l=\frac{1}{\omega-l\omega_B}\mathcal{Y}_l, \end{equation}

where

(3.35)\begin{equation} \mathcal{Y}_l(x_l,\zeta_T)=\sum_{s=-\infty}^{\infty}\int_{-\infty}^{\infty} \,{\rm d}v_\| F_\| \frac{{\rm J}_s^2(\zeta_Tv_T/v_\|)}{1-k_{\|,s}v_\|/(k_\|v_Tx_l)}, \end{equation}

$x_l =(\omega -l\bar \omega _B)/(k_\| v_T)$, $\zeta _T=lqr/\rho _T$.

Due to (3.30), the eigenmode equation in the considered case can be obtained from (3.21) by replacing $Z_l(x_l)$ with $\mathcal {Y}_l(x_l,\zeta _T)$. To draw this conclusion we substituted (3.30) into the angle-averaged (3.9) and then used (3.20) with flux surface averaged right-hand side. Thus, $\mathcal {Y}_l(x_l,\zeta _T)$ plays the role of the dispersion function $Z_l(x_l)$ when particle gyration is affected by toroidicity. Its features are considered in appendix A.

3.3. Equation for subharmonic modes

Based on (2.18), in the previous section we neglected the toroidal drift. However, finite Larmor radius and non-uniform particle gyration do not manifest themselves in the $l=0$ term of the eigenmode equations, which can describe subharmonic modes ($\omega \lesssim \omega _{Bi}$). Therefore, a question arises about the role of finite width of drift orbits for these modes. Below we consider this issue.

We proceed from (3.9) with $A_\nu (\tau, t)$ given by (3.8) which in the considered case is

(3.36)\begin{equation} A_\nu (t,\tau) = \exp\left [-{\rm i}\int_t^{\tau} \,{\rm d}t^\prime (\omega -k_\|v_\| - \boldsymbol{k}_{{\perp},\nu} \boldsymbol{v}_D)\right ], \end{equation}

where $\boldsymbol {k}_{\perp,\nu } \boldsymbol {v}_\perp =k_{\perp,\nu } v_D \sin (\vartheta +\psi )$. Using $\dot {\vartheta }=v_\|/(qR)$ and

(3.37)\begin{equation} \exp({{\rm i}\xi_D\cos\vartheta_\tau}) ={\rm i}^s\sum_{s=-\infty}^{\infty} {\rm J}_s(\xi_D) {\rm e}^{{\rm i}s\vartheta_\tau}, \end{equation}

with $\vartheta _\tau =\vartheta (\tau )+\psi$, $\xi _D=k_\perp v_DqR/v_\|$, we obtain

(3.38)\begin{equation} A_\nu (t,\tau) = \sum_s {\rm i}^s{\rm J}_s(\xi_{D,\nu})\exp\left({{\rm i}\int_t^{\tau} \,{\rm d}t^\prime \varOmega_{0,s}}\right) \exp({{\rm i}s\vartheta_t -{\rm i}\xi_{D,\nu}\cos\vartheta_t}), \end{equation}

where $\varOmega _{0,s}=\omega -k_{\|,s}v_\|$. The flux surface averaging with

(3.39)\begin{equation} {\rm J}_p(x)={\rm i}^p \oint \frac{{\rm d}\chi}{2{\rm \pi}} \exp({-{\rm i}x\cos\chi+{\rm i}p\chi}) \end{equation}

and time integration yield

(3.40)\begin{equation} \int_{-\infty}^t\,{\rm d}\tau\oint \frac{{\rm d}\vartheta}{2{\rm \pi}} A_\nu (t,\tau)={\rm i}\sum_{s=-\infty}^{\infty} \frac{{\rm J}_s^2(\xi_D)}{\varOmega_{0,s}}. \end{equation}

Using (3.40), (3.9) and (3.20) we obtain

(3.41)\begin{align} & \frac{1}{r}\frac{{\rm d}}{{\rm d}r}r\frac{{\rm d} \hat{\varPhi}}{{\rm d} r} - \left (\frac{m^2}{r^2} + \frac{n^2}{R^2}\right )\hat{\varPhi} \nonumber\\ & \qquad =\sum_{j=e,i}\sum_{\nu=-\infty}^{\infty}\frac{1}{d_j^2}\left \langle 1 - \sum_{s=-\infty}^{\infty}{\rm J}_s^2(\xi_{D,\nu})\frac{\omega}{\varOmega_{0,s}} \right \rangle\varPhi_\nu {\rm e}^{{\rm i}k_\nu r}, \end{align}

where $\langle (\cdots )\rangle \equiv \int \,{\rm d}^3v (\cdots )F_\|F_\perp$.

The mode frequency below but close to $\omega _{Bi}$ well exceeds $v_\| /(qR)$, by a factor of $qR/\rho$. On the other hand, when $\xi _D \sim 1$ or less, a few terms in the sum $\sum _{-\infty }^{\infty } {\rm J}_s^2(\xi _D)$ are sufficient to make this sum very close to unity, $\sum _{-s_m}^{s_m} {\rm J}_s^2(\xi _{D,\nu })\approx 1$. Due to this, the right-hand side of (3.41) reduces to

(3.42)\begin{equation} \sum_{j=e,i}\frac{1}{d_j^2}[1-Z_0(x_{0,j})] \hat\varPhi . \end{equation}

Thus, the influence of the toroidal drift on the equation of modes with $\omega \lesssim \omega _{Bi}$ is negligible.