On equations for ion cyclotron modes in ‘warm’ bounded plasmas

Abstract

An equation describing eigenmodes in the ion cyclotron frequency range in ‘warm’ bounded plasmas, i.e. eigenmodes which are absent in the two-fluid model but exist in kinetic theory due to finite Larmor radius of the ions, is derived for the first time. It is valid for electrostatic modes but the developed approach is generic. Calculations are carried out for two cases: first, for a homogeneous magnetic field; second, taking into account the effects of toroidicity in tokamaks. It is found that, in general, equations for eigenmodes in warm plasmas do not reduce to second-order differential equations (in contrast to those which are usually used to describe the radial structure of eigenmodes in fusion devices). The study of modes in warm plasmas is of interest, in particular, in connection with the recent observations of superthermal ion cyclotron emission in the NSTX-U spherical torus and DIII-D tokamak, which can be hardly explained by conventional theories employing fast magnetoacoustic modes.

1. Introduction

Eigenmodes in the ion cyclotron frequency range in toroidal plasmas are usually described within a cold-plasma approximation (two-fluid model), in which case only fast magnetoacoustic modes (FMM) (known also as CAE (compressional Alfvén eigenmodes)) exist. However, recent experiments with ion cyclotron emission (ICE) on the NSTX-U spherical torus (Fredrickson et al. 2021) and DIII-D tokamak (Crocker et al. 2022) indicate that this approximation can be unsuitable. First, it was found that the emission frequency in NSTX-U does not follow an Alfvénic scaling with density; this newly discovered type of chirping ICE does not follow the well-known ICE phenomenology on JET, TFTR, KSTAR or LHD. Second, a new diagnostic with Doppler backscattering on the DIII-D has shown that the wavenumber exceeded that of FMM: $k\gg \omega /v_A$, where $k$ and $\omega$ are the wavenumber and wave frequency, respectively, and $v_A$ the Alfvén velocity (Crocker et al. 2022). Therefore, it was concluded in Crocker et al. (2022) that the observed ICE resulted from destabilized ion cyclotron waves, either electrostatic waves (which presumably took place in TFTR (Dendy et al. 1994)) or electromagnetic waves. In contrast to FMM and Alfvén waves, these modes are absent in ideal magnetohydrodynamics (MHD) and the two-fluid model with zero temperature. They require a kinetic theory for their description in order to take into account finite Larmor radius of thermal ions of the bulk plasma. We therefore refer to them as modes in ‘warm’ plasmas.

The results of theories of waves in warm plasmas can be found in, for example, overviews (Akhiezer et al. 1975; Stix 1992) and references therein. However, they are obtained for infinite plasmas and for this reason the developed theories cannot determine the eigenmode location and its radial structure in warm plasmas; they can be useful for a local analysis of possible instabilities, but are not sufficient for a reliable identification of destabilized modes. On the other hand, conventional equations for eigenmodes in bounded plasmas (the second-order linear differential equations) can be inappropriate. The reason is that coefficients in the eigenmode equations in warm plasmas depend on the mode structure. Moreover, eigenmodes are standing waves, whereas the known way to take into account finite Larmor radius employs travelling waves.

The mentioned experimental findings, in particular $k \gg \omega /v_A$, are consistent with the high-$k$ part of the FMM branch of the dispersion equation. However, this does not mean that other modes could not be responsible for these findings. In particular, electrostatic modes, the equations for which will be derived in this work, well satisfy $k \gg \omega /v_A$. Although the radial mode structure was determined experimentally, no theoretical calculations of radial dependence of the mode amplitudes of cyclotron waves were carried out. Accordingly, the experimentally observed modes were not conclusively identified.

All these facts motivated us to carry out this work devoted to ICE eigenmode equations in a warm plasma.

The work is organized as follows. Section 2 outlines our approach to the description of eigenmodes. It consists of two parts. The first part (§ 2.1) deals with the peculiarities of equations for eigenmodes in bounded warm plasmas and suggests a way to obtain equations for cyclotron modes. The second part (§ 2.2) analyses key parameters associated with toroidicity and discusses their roles. In § 3 an equation for electrostatic eigenmodes valid for the cyclotron frequency range is derived: in § 3.1 the approximation of a homogeneous magnetic field is used; an effect of toroidicity – non-uniform cyclotron rotation (gyration) of the particles – is taken into account in § 3.2. In addition, an equation for subharmonic modes in tokamaks is derived in § 3.3. Section 4 summarizes the results obtained. Appendix A considers features of the dispersion function associated with non-uniform particle gyration.

2. Approach to the description of eigenmodes

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 (1992)):

(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 (1993) and the overview Gorelenkov (2016). 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 (1993). 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 1993),

(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 (1993), we consider the JET preliminary tritium experiment where the ICE was observed (Cottrell et al. 1993), 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 (1993) 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 2020). 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. 1997) observed in JET (Cottrell et al. 1993)).

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. Derivation of equations for electrostatic eigenmodes

In this section we derive an eigenmode equation assuming $\delta \boldsymbol {B} =0$, in which case a scalar potential ($\delta \varPhi$) describes the perturbed electric field, $\delta \boldsymbol {E} =-\boldsymbol {\nabla }\delta \varPhi$. In order to calculate the disturbed distribution function of the particles we will expand $\delta \varPhi$ in a Fourier series, as outlined in § 2.1.

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. 2017). 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 1982; Mikhailovskii 1986). 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 1967). Note that $Z_l(x) =-xZ_{dis}(x)$, where $Z_{dis}(x)$ is the plasma dispersion function of Richardson (2019).

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 (1961) 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 1986; Belikov & Kolesnichenko 1982)

(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.

4. Summary

Cyclotron waves can manifest themselves in fusion devices and in space plasmas (Dendy 1994). However, minor attention was paid to eigenmodes in warm plasmas after the first works carried out in the middle of the last century. In particular, no attempts to describe the radial structure of those eigenmodes were done yet. This can be explained by the fact that there are many experimental observations of destabilized MHD modes, especially toroidicity-induced Alfvén eigenmodes, and at the same time the theory of modes in cold plasma is simpler (because they are described by MHD and two-fluid models and do not require kinetic theory). This work represents an indispensable step towards a study of eigenmodes in warm bounded plasmas. It was motivated, first of all, by the ICE experiments on NSTX-U and DIII-D (Fredrickson et al. 2021; Crocker et al. 2022), which can be hardly described by existing theories employing FMMs. Some other recent DIII-D results on ICE can be found in Thome et al. (2019), DeGrandchamp et al. (2021) and DeGrandchamp et al. (2022).

The results obtained within the work can be summarized as follows.

A method for derivation of equations for the mentioned eigenmodes, based on expansion of the radial profile of the mode amplitude in a Fourier series with coefficients in the form of plane waves is proposed. Using this method, an equation for electrostatic modes with frequencies in the range of ion cyclotron harmonics in a Maxwellian plasma is derived for the first time, see (3.21). In fact, this is an integro-differential equation. It reduces to a second-order differential equation in the special case of $k_r\ll k$.

Equation (3.21) is derived in the assumption that the magnetic field is homogeneous. A comprehensive treatment of toroidicity would have resulted in an equation which has a form hardly suitable for applications (more details are given in § 2.2). Therefore, we restricted ourselves to taking into account the non-uniform gyration of particles, which seems to represent the main effect of toroidicity on eigenmodes when toroidicity is small ($\sqrt {\epsilon } \ll 1$) and/or $k_\|qR\ll 1$. It is found that non-uniform gyration leads to replacement of the known dispersion function $Z_l(x_l)$ in the mode equation with the function $\mathcal {Y}_l(x_l,\zeta _T)$. Both these functions are proportional to $|\omega -l\bar {\omega }_B|^2$ when $|\omega -l\bar {\omega }_B|$ is sufficiently small, and weakly depend on $|\omega -l\bar {\omega }_B|$ in the opposite case.

An equation for subharmonic modes ($\omega \lesssim \omega _{Bi}$), that takes into account finite drift-orbit width, is derived. It is found that toroidicity is not important for these modes.

A by-product of this work is finding a minimum magnitude of the poloidal mode number in the known ICE model (Coppi 1993), which restricts applicability of that model to high cyclotron harmonics.