2.1 Slab geometry

Following Helander (Reference Helander2014) and Helander & Connor (Reference Helander and Connor2016), it is convenient to write the gyrokinetic distribution function in the form

(2.1)\begin{equation} f_a = f_{a0} \left(1 - \frac{e_a \langle \phi \rangle}{T_a}\right) + g_a = f_{a0} + f_{a1},\quad f_{a1} ={-}\frac{e_a \langle \phi \rangle}{T_a}f_{a0} + g_a . \end{equation}

Here, $f_{a0}$ is a Maxwellian, $a$ is the species index with $a=e$ corresponding to electrons, $a = p$ to positrons and $a = i$ to the ions, $e_a$ is the electric charge, $f_{a1}$ is the perturbed part of the distribution function and $g_a$ is the non-adiabatic part of $f_{a1}$. The linearized gyrokinetic equation in this notation is

(2.2)\begin{equation} {\rm i} v_{\|} \boldsymbol{\nabla}_{\|} g_a + (\omega - \omega_{{\rm d}a}) g_a = \frac{e_a}{T_a} {\rm J}_0\left(\frac{k_{{\perp}}v_{{\perp}}}{\omega_{{\rm c}a}}\right) \left(\omega - \omega^{{\rm T}}_{*a}\right) (\phi - v_{\|} A_{\|}) f_{a0} \end{equation}

with ${\rm J}_0$ the Bessel function, $\omega$ the complex frequency of the mode, $\omega _{{\rm c}a}$ the cyclotron frequency, $k_{\perp }$ the component of the wavenumber perpendicular to the ambient magnetic field, $v_{\|}$ and $v_{\perp }$ the parallel and perpendicular velocities, $\phi$ the perturbed electrostatic potential and $A_{\|}$ the perturbed parallel magnetic potential in the Coulomb gauge. The plasma pressure has been assumed small enough that magnetic-field fluctuations parallel to the equilibrium field can be neglected. We consider an unsheared slab geometry with coordinates $(x, y, z)$, a uniform magnetic field ${\boldsymbol {B}} = B {\boldsymbol {e}}_z$ pointing in the $z$-direction and plasma profiles which are non-uniform in the $x$-direction. In the slab geometry, the drift frequency $\omega _{{\rm d}a} = 0$. Other notations used are

(2.3a–c)\begin{gather} \omega_{*a}^{{\rm T}} = \omega_{*a} \left[ 1 + \eta_a \left( \frac{v^2}{v_{{\rm th}a}^2} - \frac{3}{2}\right) \right],\quad v = \sqrt{v_{\|}^2 + v_{{\perp}}^2},\quad k_{{\perp}} = \sqrt{k_x^2 + k_y^2} \end{gather}

(2.4a–d)\begin{gather}\omega_{*a} = \frac{k_y T_a}{e_a B} \frac{{\rm d}{\ln n_a}}{{\rm d}{x}},\quad \eta_a = \frac{{\rm d}{\ln T_a}}{{\rm d}{\ln n_a}},\quad v_{{\rm th}a} = \sqrt{\frac{2 T_a}{m_a}},\quad \omega_{{\rm c}a} = \frac{e_a B}{m_a}. \end{gather}

Here, $m_a$ is the particle mass, $n_a$ is the ambient particle density and the sign convention is such that $\omega _{*i} \le 0$, $\omega _{*p} \le 0$, and $\omega _{*e} \ge 0$ for ${\rm d} \ln n_a / {\rm d} x \le 0$. For simplicity, we will assume $k_x = 0$ and $k_{\perp } = k_y$ throughout the paper. Taking the Fourier transform along the parallel coordinate, we obtain

(2.5)\begin{equation} ( \omega - k_{\|} v_{\|} ) g_a = \frac{e_a}{T_a} {\rm J}_0\left(\frac{k_{{\perp}}v_{{\perp}}}{\omega_{{\rm c}a}}\right) \left(\omega - \omega^{{\rm T}}_{*a}\right) (\phi - v_{\|} A_{\|}) f_{a0}. \end{equation}

This equation is trivially solved

(2.6)\begin{equation} g_a = \frac{\omega - \omega^{{\rm T}}_{*a}}{\omega - k_{\|} v_{\|}} \frac{e_a f_{a0}}{T_a} {\rm J}_0 (\phi - v_{\|}A_{\|}). \end{equation}

The gyrokinetic quasineutrality condition and the parallel Ampere's law are

(2.7a,b)\begin{equation} \left( \sum_a \frac{n_a e_a^2}{T_a} + \epsilon_0 k_{{\perp}}^2 \right) \phi = \sum_a e_a \int g_a \,{\rm J}_0 \,{\rm d}^3 v,\quad A_{\|} = \frac{\mu_0}{k_{{\perp}}^2} \sum_a e_a \int v_{\|} g_a \,{\rm J}_0 \,{\rm d}^3 v. \end{equation}

Here, $\epsilon _0$ is the electric permittivity and $\mu _0$ is the magnetic permeability of vacuum. For the electromagnetic dispersion relation, it is convenient to define

(2.8)\begin{equation} W_{n a} ={-}\frac{1}{n_a v_{{\rm th}a}^n} \int \frac{\omega - \omega^{{\rm T}}_{*a}}{\omega - k_{\|} v_{\|}}\,{\rm J}_0^2 f_{a0} v_{\|}^n \, {\rm d}^3 v. \end{equation}

Taking velocity-space integrals, one finds

(2.9)\begin{equation} W_{na} = \zeta_a \left\{ \left( 1 - \frac{\omega_{*a}}{\omega} \right) Z_{na} \varGamma_{0a} + \frac{\omega_{*a} \eta_a}{\omega} \left[\frac{3}{2} Z_{na} \varGamma_{0a} - Z_{na} \varGamma_{*a} - Z_{n+2,a} \varGamma_{0a} \right]\right\}. \end{equation}

Here, the following notation is employed:

(2.10a–d)\begin{gather} \frac{1}{\lambda_{{\rm D}a}^2} = \frac{e_a^2 n_a}{\epsilon_0 T_a},\quad \frac{1}{\lambda_{\rm D}^2} = \sum_a \frac{1}{\lambda_{{\rm D}a}^2},\quad b_a = k_{{\perp}}^2 \rho_a^2,\quad \rho_a = \frac{\sqrt{m_a T_a}}{|e_a| B} \end{gather}

(2.11)\begin{gather}\varGamma_{*a} = \varGamma_{0a} - b_a \left[\varGamma_{0a} - \varGamma_{1a}\right],\quad \varGamma_{0a} = I_0(b_a) \,{\rm e}^{{-}b_a},\quad \varGamma_{1a} = I_1(b_a) \,{\rm e}^{{-}b_a} \end{gather}

(2.12)\begin{gather}Z_{na} = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^{\infty} \frac{x^n \,{\rm e}^{{-}x^2} \,{\rm d} x}{x - \zeta_a},\quad \zeta_a = \frac{\omega}{k_{\|} v_{{\rm th}a}}, \end{gather}

with $I_0$ and $I_1$ denoting the modified Bessel functions of the first kind. Note the presence of the Debye scale in (2.7a,b) and (2.10a–d) usually ignored in gyrokinetic applications for fusion plasmas already at the level of the phase-space Lagrangian. The resulting electromagnetic dispersion relation in slab geometry reads

(2.13)\begin{align} & \left( 1 + k_{{\perp}}^2 \lambda_D^2 + \sum_a \frac{\lambda_D^2}{\lambda_{Da}^2} W_{0a} \right) \left( 1 - 2 \sum_a \frac{\beta_a}{k_{{\perp}}^2 \rho_a^2} W_{2a} \right) \nonumber\\ & \quad +2 \sum_a \frac{\lambda_D^2}{\lambda_{Da}^2} W_{1a} v_{{\rm th}a} \sum_a\frac{\beta_a}{k_{{\perp}}^2 \rho_a^2} \frac{W_{1a}}{v_{{\rm th}a}} = 0. \end{align}

Here, $\beta _a = \mu _0 n_a T_a / B^2$. The electrostatic limit corresponds, as usual, to $\beta _a = 0$.

This dispersion relation has been solved in Mishchenko et al. (Reference Mishchenko, Zocco, Helander and Könies2018b). It was found that pair plasmas can support the gyrokinetic ion-temperature-gradient-driven (ITG), electron-temperature-gradient-driven (ETG) and universal instabilities even in a slab geometry if the proton fraction exceeds some threshold. In practice, however, this threshold is usually quite large, hopefully large enough to keep the proton content below this value in pair-plasma experiments (Pedersen et al. Reference Pedersen, Danielson, Hugenschmidt, Marx, Sarasola, Schauer, Schweikhard, Surko and Winkler2012). These results extend the finding of Helander (Reference Helander2014) that pair plasmas are stable to gyrokinetic modes in the absence of magnetic curvature for the cases with small to moderate proton contamination. We find, however, that pure pair plasmas can have temperature-gradient-driven instabilities, if the electron and the positron temperature profiles differ. In reality, such profiles are unlikely in steady state, since the characteristic time of energy exchange between the species is comparable to the equilibration time of each species. Generalization of the local dispersion relation to the case of non-neutral plasma is straightforward (Kennedy & Mishchenko Reference Kennedy and Mishchenko2019) providing that effects of a strong electric field, normally existing in non-neutral plasmas, can be cast in the form of a simple Doppler shift for the frequency. The dispersion relation for this shifted frequency coincides with (2.13).

2.2 Dipole geometry

In a dipole geometry, the drift frequency $\omega _{{\rm d}a} \ne 0$, so that the gyrokinetic equation (2.2) has to be used where now

(2.14)\begin{gather} \omega_{{\rm d}a} = {\boldsymbol{k}}_{{\perp}} \boldsymbol{\cdot} {\boldsymbol{v}}_{{\rm d}a},\quad {\boldsymbol{v}}_{{\rm d}a} = \left(m v_{\|}^2 + \mu B \right) \frac{{\boldsymbol{b}} \times \boldsymbol{\nabla} B}{e_a B^2},\quad \mu = \frac{m_a v_{{\perp}}^2}{2 B} \end{gather}

(2.15a,b)\begin{gather}{\boldsymbol{k}}_{{\perp}} = k_{\psi} \boldsymbol{\nabla} \psi + k_{\varphi} \boldsymbol{\nabla} \varphi,\quad \omega_{*a} = \frac{k_{\varphi} T_a}{e_a} \frac{{\rm d}{\ln n_a}}{{\rm d}{\psi}} , \end{gather}

with $\psi$ the poloidal flux and $\varphi$ the polar (toroidal) angle. We will assume electrostatic perturbations ($A_{\|} = 0$) and the drift-kinetic limit, i.e. $k_{\perp } v_{{\rm th}a}/\omega _{{\rm c}a} \ll 1$ so that ${\rm J}_0 \approx 1$.

Expanding the distribution function $g_a = g^{(0)}_a + g^{(1)}_a + \cdots$ in the small parameter $\varepsilon _b = \omega /\omega _b$ with $\omega _b$ the bounce frequency, we obtain in the lowest order

(2.16)\begin{equation} v_{\|} \boldsymbol{\nabla}_{\|} g_a^{(0)} = 0, \end{equation}

implying that $g_a^{(0)}$ coincides with its bounce average, $g_a^{(0)} = \bar {g}_a^{(0)}$, where the bounce-average operation is defined as

(2.17)\begin{equation} \overline{(\cdots)} = \oint (\cdots) \frac{{\rm d} l}{v_{\|}}/\oint \frac{{\rm d} l}{v_{\|}}. \end{equation}

Here, $l$ is the arc length measured along a magnetic-field line and the integration is performed between bounce points for trapped particles, and over the entire closed field line for passing particles. Applying the bounce average at the next order in $\varepsilon _b$, we obtain

(2.18)\begin{equation} (\omega - \bar{\omega}_{{\rm d}a}) g^{(0)}_a = (\omega - \omega_{{\ast} a}^{{\rm T}}) \frac{e_a \bar{\phi}}{T_a} f_{a0}, \end{equation}

where we have neglected magnetic fluctuations. We assume the background temperature and the density profiles of the electrons and the positrons to be identical, and invoke the Poisson equation for the charge density perturbations

(2.19)\begin{equation} \left( \sum_{a=e,p} \frac{n_a e_a^2}{T_a} + \epsilon_0 k_{{\perp}}^2 \right) \phi = \sum_{a=e,p} e_a \int g^{(0)}_a \,{\rm d}^3 v . \end{equation}

We find that the perturbed electrostatic potential satisfies the equation

(2.20)\begin{equation} \left(1 + k_{{\perp}}^2 \lambda_D^2 \right) \phi = \frac{1}{n_0} \int \frac{\omega^2 - \bar{\omega}_d \omega_*^{{\rm T}}}{\omega^2 - \bar{\omega}_d^2} \bar{\phi} f_0 \,{\rm d}^3 v. \end{equation}

Here, and in the following, we use the notation $\omega _*^{{\rm T}} \equiv \omega _{*e}^{{\rm T}}$, $\omega _* \equiv \omega _{*e}$, $\bar {\omega }_d \equiv \bar {\omega }_{{\rm d}e}$, $n_0 = n_e$, $T_0 \equiv T_e$ and the Debye length is defined as usual, $\lambda _D = \sqrt {\epsilon _0 T_0 / (2 n_0 \,{\rm e}^2)}$.

We have solved this equation in Mishchenko et al. (Reference Mishchenko, Plunk and Helander2018a) where the drift-kinetic stability of a pair plasma confined by a dipole magnetic field has been studied. It has been found that pair plasmas can be unstable in a dipole geometry even for perfectly coinciding electron and positron profiles, in the absence of any contamination, and for quasineutral plasmas. The reason for the instability is related to the fact that the curvature drift depends on the sign of the particle charge, which is opposite in the case of the electrons and positrons. A detailed study of instabilities in dipole pair plasmas has been carried out in Mishchenko et al. (Reference Mishchenko, Plunk and Helander2018a). In contrast, one needs some violation of the symmetry between the species for instability in a slab.

3.1 Relativistic gyrokinetic equation

The relativistic gyrokinetic equation has been derived in Brizard & Chan (Reference Brizard and Chan1999) as

(3.1)\begin{equation} \frac{\partial{f_a}}{\partial{t}} + \left(\frac{{\boldsymbol{B}}^*}{B_{\|}^*} \frac{\partial{H_{gy}}}{\partial{p_{\|}}} + c \frac{{\boldsymbol{b}}\times\boldsymbol{\nabla} H_{gy}}{e_a B_{\|}^*} \right) \boldsymbol{\cdot}\boldsymbol{\nabla} f_a - \frac{{\boldsymbol{B}}^*\boldsymbol{\cdot}\boldsymbol{\nabla} H_{gy}}{B_{\|}^*} \frac{\partial{f_a}}{\partial{p_{\|}}} = 0. \end{equation}

Note that the CGS unit system is employed throughout this section since it is more convenient than the SI (MKS) units in relativistic calculations. To the first order in the perturbation amplitude, the gyrokinetic relativistic Hamiltonian is given by the expression

(3.2)\begin{equation} H_{gy} = \gamma m_a c^2 + e_a \left\langle \phi - \frac{p_{\|}}{\gamma m_a c} A_{\|} - \frac{1}{\gamma c} \sqrt{\frac{2 \mu B}{m_a}} {\pmb{\perp}}\boldsymbol{\cdot}{\boldsymbol{A}} \right\rangle, \end{equation}

with ${\pmb {\perp }}$ the unit vector directed along the gyro-motion of the particle and the gyro-average operation $\langle \phi \rangle = \oint \phi ({\boldsymbol {R}}+{\boldsymbol {\rho }}) \,{\rm d}\theta /(2{\rm \pi} )$ defined as usual for the gyro-radius

(3.3)\begin{equation} {\boldsymbol{\rho}} = \frac{1}{\omega_{Ba}} \sqrt{\frac{2 \mu B}{m_a}} {\boldsymbol{\zeta}},\quad \omega_{Ba} = \frac{e_a B}{m_a c}, \end{equation}

where ${\boldsymbol {\zeta }} = ({\boldsymbol {e}}_1 \cos \theta - {\boldsymbol {e}}_2 \sin \theta )$ is the unit vector directed along the gyro-radius rotating in the fixed basis $({\boldsymbol {e}}_1, {\boldsymbol {e_2}}, {\boldsymbol {b}})$ and $\omega _{Ba}$ is the rest-mass gyro-frequency. Following Brizard & Chan (Reference Brizard and Chan1999), the gyrokinetic relativistic Lorentz factor is

(3.4)\begin{equation} \gamma = \sqrt{1 + \frac{p^2}{m^2 c^2}} = \sqrt{1 + \frac{2 \mu B}{m_a c^2} + \left(\frac{p_{\|}}{m_a c}\right)^2},\quad p^2 = p_{\|}^2 + 2 m_a \mu B, \end{equation}

and

(3.5a–c)\begin{equation} {\boldsymbol{B}}^* = {\boldsymbol{B}} + c \frac{p_{\|}}{e_a} \boldsymbol{\nabla}\times{\boldsymbol{b}},\quad B_{\|}^* = {\boldsymbol{b}}\boldsymbol{\cdot}{\boldsymbol{B}}^*,\quad {\boldsymbol{b}}^* = \frac{{\boldsymbol{B}}^*}{B_{\|}^*}. \end{equation}

We can write

(3.6)\begin{gather} \frac{\partial{H_{gy}}}{\partial{p_{\|}}} = \left(\frac{p_{\|}}{m_a \gamma} - \frac{e_a}{m_a}\frac{\langle A_{\|} \rangle}{\gamma c}\right) + \frac{e_a p_{\|}}{m_a^3 c^3 \gamma^3} \left[p_{\|}\langle A_{\|} \rangle + \sqrt{2 m_a \mu B} \left\langle {\pmb{\perp}}\boldsymbol{\cdot}{\boldsymbol{A}} \right\rangle\right] \end{gather}

(3.7)\begin{align} \boldsymbol{\nabla} H_{gy} & = \left[m c^2 + \frac{e_a}{\gamma^2 m_a c}\left(p_{\|}\langle A_{\|} \rangle + \sqrt{2 m_a \mu B} \left\langle {\pmb{\perp}}\boldsymbol{\cdot}{\boldsymbol{A}} \right\rangle\right)\right] \boldsymbol{\nabla}\gamma - \frac{e_a}{\gamma c} \sqrt{\frac{2 \mu B}{m_a}} \frac{\boldsymbol{\nabla} B}{2 B} \left\langle {\pmb{\perp}}\boldsymbol{\cdot}{\boldsymbol{A}} \right\rangle \nonumber\\ & \quad + e_a \left(\boldsymbol{\nabla}\langle \phi \rangle - \frac{p_{\|}}{\gamma m c} \boldsymbol{\nabla}\langle A_{\|} \rangle - \frac{\sqrt{2 m_a \mu B}}{\gamma m c}\boldsymbol{\nabla}\left\langle {\pmb{\perp}}\boldsymbol{\cdot}{\boldsymbol{A}} \right\rangle\right),\quad \boldsymbol{\nabla}\gamma = \frac{\mu\boldsymbol{\nabla} B}{m c^2 \gamma}. \end{align}

3.2 Electrostatic dispersion relation

For electrostatic waves ${\boldsymbol {A}} = 0$. In this case, the linearized relativistic gyrokinetic equation in a slab geometry takes the form

(3.8)\begin{equation} \frac{\partial{f_{1a}}}{\partial{t}} + \frac{p_{\|}}{m_a \gamma} {\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} f_{1a} = e_a {\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\langle \phi \rangle \frac{\partial{F_{0a}}}{\partial{p_{\|}}} - c \frac{{\boldsymbol{b}}\times\boldsymbol{\nabla}\langle \phi \rangle}{B}\boldsymbol{\cdot}\boldsymbol{\nabla} F_{0a}. \end{equation}

For the background, we can choose the Maxwell–Jüttner distribution function describing the relativistic ideal gas, see Jüttner (Reference Jüttner1911), Zenitani (Reference Zenitani2015) and Cercignani & Kremer (Reference Cercignani and Kremer2012),

(3.9)\begin{equation} F_{0a} = N_a \exp\left[-\frac{c}{T_a} \sqrt{p^2 + m_a^2 c^2}\right],\quad N_a = \frac{n_{0a}}{4{\rm \pi} m_a^3 c^3} \left[\frac{T_a}{m_a c^2} K_2\left(\frac{m_a c^2}{T_a}\right)\right]^{{-}1}. \end{equation}

Here, $K_2(z)$ is the modified Bessel function of the second kind which has the asymptotic behaviour

(3.10)\begin{equation} K_2(z) \approx \sqrt{\frac{\rm \pi}{2 z}} {\rm e}^{{-}z},\quad z \gg1. \end{equation}

Using this asymptotic expression, we can recover the usual Maxwellian distribution function for $T_a \ll m c^2$ and $p \ll m c$. The distribution function (3.9) neglects the pair formation and quantum effects but it can be used as an approximation for already created relativistic electron–positron plasmas. Taking derivatives, we obtain

(3.11a,b)\begin{gather} \frac{\partial{F_{0a}}}{\partial{p_{\|}}} ={-}\frac{p_{\|}}{m_a \gamma} \frac{F_{0a}}{T_a},\quad \frac{\partial{F_{0a}}}{\partial{\mu}} ={-}\frac{B}{\gamma} \frac{F_{0a}}{T_a} \end{gather}

(3.12)\begin{gather}\boldsymbol{\nabla} F_{0a} = \left[\frac{\boldsymbol{\nabla} N_a}{N_a} + \frac{\gamma m_a c^2}{T_a} \left(\frac{\boldsymbol{\nabla} T_a}{T_a} - \frac{\boldsymbol{\nabla} \gamma}{\gamma}\right) \right] F_{0a}. \end{gather}

Note that $\boldsymbol {\nabla } \gamma$ will contribute only the case of an inhomogeneous magnetic field. For a uniform magnetic field considered in this section, $\boldsymbol {\nabla } \gamma = 0$. Using these expressions, the relativistic gyrokinetic equation can be written in the form

(3.13)\begin{equation} \frac{\partial{g_a}}{\partial{t}} + \frac{p_{\|}}{m_a \gamma} \boldsymbol{\nabla}_{\|} g_a = \frac{e_a F_{0a}}{T_a} \left( \frac{\partial{\langle \phi \rangle}}{\partial{t}} - \frac{c T_a}{e_a B} \frac{\boldsymbol{\nabla} F_{0a} \times {\boldsymbol{b}}}{F_{0a}} \boldsymbol{\cdot}\boldsymbol{\nabla}\langle \phi \rangle \right), \end{equation}

with the usual definition of the non-adiabatic part of the distribution function $g_a$

(3.14)\begin{equation} f_a = F_{0a} \left(1 - \frac{e_a \langle \phi \rangle}{T_a}\right) + g_a = F_{0a} + f_{1a},\quad f_{1a} ={-}\frac{e_a \langle \phi \rangle}{T_a}F_{0a} + g_a. \end{equation}

For plane waves, $g_a = \hat {g}_{a} \exp ({\rm i} \omega t - {\rm i} {\boldsymbol {k}}\boldsymbol {\cdot }{\boldsymbol {x}})$ and $\phi = \hat {\phi } \exp ({\rm i} \omega t - {\rm i} {\boldsymbol {k}}\boldsymbol {\cdot }{\boldsymbol {x}})$, we can write

(3.15)\begin{equation} \left(\omega - \frac{k_{\|} p_{\|}}{m_a\gamma}\right) \hat{g}_a = \left(\omega - \omega_{*a}^{{\rm T}}\right) \frac{e_a \hat{\phi}}{T_a} F_{0a} {\rm J}_0(k_{{\perp}}\rho),\quad \omega_{*a}^{{\rm T}} = \frac{c T_a}{e_a B}({\boldsymbol{k}} \times {\boldsymbol{b}}) \boldsymbol{\cdot} \frac{\boldsymbol{\nabla} F_{0a}}{F_{0a}}. \end{equation}

Following Brizard & Chan (Reference Brizard and Chan1999), we write the gyrokinetic Maxwell equations in the low-frequency limit neglecting the displacement current

(3.16)\begin{gather} -\frac{1}{4{\rm \pi}} \nabla_{{\perp}}^2 A^{\alpha} = \sum_a e_a \int {\rm d}^6 Z \delta^3_{gy} \frac{p^{\alpha}}{m_a \gamma c} \left[f_{1a} + \{S_1, F_{0a}\} + \frac{e_a {\boldsymbol{A}}}{c}\{{\boldsymbol{R}} + {\boldsymbol{\rho}}, F_{0a}\}\right] \end{gather}

(3.17a–c)\begin{gather}S_1 = \frac{e_a \gamma}{\omega_{Ba}} \int^{\theta} \tilde{\psi} \, {\rm d}\theta',\quad \psi = \phi - \frac{{\boldsymbol{p}}}{\gamma m c} \boldsymbol{\cdot} {\boldsymbol{A}},\quad \{S_1, F_{0a}\} \approx \frac{e_a}{m_a c} \frac{\partial{S_1}}{\partial{\theta}} \frac{\partial{F_{0a}}}{\partial{\mu}} \end{gather}

(3.18)\begin{gather}\tilde{\psi} = \psi({\boldsymbol{R}}+{\boldsymbol{\rho}}) - \langle \psi \rangle,\quad \langle \psi \rangle = \oint \psi({\boldsymbol{R}}+{\boldsymbol{\rho}}) \frac{{\rm d}\theta}{2{\rm \pi}} , \end{gather}

with $A^{\alpha } = (\phi, {\boldsymbol {A}})$ the perturbed four-potential, $p^{\alpha } = (\gamma m_a c, {\boldsymbol {p}})$ the four-momentum, ${\boldsymbol {p}} = p_{\|} {\boldsymbol {b}} + \sqrt {2 m_a \mu B} {\pmb {\perp }}$, ${\rm d}^6 Z = m_a B_{\|}^* \,{\rm d}{\boldsymbol {R}} \,{\rm d} p_{\|} \,{\rm d}\mu \,{\rm d}\theta$ and $\delta ^3_{gy} = \delta ({\boldsymbol {R}} + {\boldsymbol {\rho }} - {\boldsymbol {x}})$.

For electrostatic perturbations, only the quasineutrality equation has to be solved

(3.19)\begin{equation} -\frac{1}{4{\rm \pi}} \nabla_{{\perp}}^2\phi = \sum_a e_a \bar{n}_{1a} - \sum_a e_a \int {\rm d}^6 Z \delta^3_{gy} \frac{e_a \tilde{\phi}}{T_a} F_{0a},\quad \tilde{\phi} = \phi - \langle \phi \rangle \end{equation}

with the perturbed density of the gyro-centres $\bar {n}_{1a} = \int {\rm d}^6 Z \,\delta ^3_{gy} f_{1a}$. In terms of the non-adiabatic part of the distribution function, the quasineutrality condition reads

(3.20)\begin{equation} -\frac{1}{4{\rm \pi}} \nabla_{{\perp}}^2\phi + \sum_a e_a \int {\rm d}^6 Z \,\delta^3_{gy} \frac{e_a \phi}{T_a} F_{0a} = \sum_a e_a \int {\rm d}^6 Z \,\delta^3_{gy} g_a. \end{equation}

For plane waves, we can write as usual

(3.21)\begin{equation} \frac{1}{4{\rm \pi}} k_{{\perp}}^2 \hat{\phi} + \sum_a \frac{e_a^2 \hat{\phi}}{T_a} n_{0a} = \sum_a e_a \int {\rm J}_0 g_a \,{\rm d}^3 p, \end{equation}

with ${\rm d}^3 p = 2{\rm \pi} m_a B_{\|}^* \,{\rm d} p_{\|} \,{\rm d}\mu$ and $n_{0a} = \int {\rm d}^3 p \,F_{0a}$. This results in the relativistic dispersion relation for gyrokinetic electrostatic waves

(3.22)\begin{equation} 1 + k_{{\perp}}^2 \lambda_D^2 + \sum_a \frac{\lambda_D^2}{\lambda_{Da}^2} W_{0a} = 0, \end{equation}

with the notation

(3.23a–c)\begin{equation} W_{0a} ={-}\frac{1}{n_{0a}} \int \frac{\omega - \omega_{*a}^{{\rm T}}}{\omega - k_{\|} p_{\|}/(m_a \gamma)} {\rm J}_0^2 F_{0a} \,{\rm d}^3 p,\quad \lambda_{Da}^2 = \frac{T_a}{4{\rm \pi} e_a^2 n_{0a}},\quad \frac{1}{\lambda_D^2} = \sum_a \frac{1}{\lambda_{Da}^2}. \end{equation}

One can easily see that the effect of spatial non-uniformity is proportional to $\omega _{*a}^{{\rm T}}$. It disappears in a pure relativistic pair plasma, similarly to the non-relativistic case, providing the electron and positron profiles are the same and magnetic drifts are absent. One can see that the resonance structure $\omega = k_{\|} p_{\|}/(m_a \gamma )$ appearing in the function $W_{0a}$ is more complex than in the non-relativistic case (recall that $\gamma$ depends both on $p_{\|}$ and the magnetic moment $\mu$). Such integrals involving the Maxwell–Jüttner distribution cannot simply be expressed through the plasma dispersion function, as they were in the non-relativistic case. Note that the relativistic generalization of the plasma dispersion functions has extensively been studied in the weakly relativistic limit (Robinson Reference Robinson1987; Castejón & Pavlov Reference Castejón and Pavlov2006) for electron cyclotron heating applications.

3.3 Electromagnetic dispersion relation

Now let us address relativistic gyrokinetic theory in the electromagnetic regime. For simplicity, we consider slab geometry and neglect the compressional component of the magnetic field perturbation assuming ${\boldsymbol {A}} \approx A_{\|}{\boldsymbol {b}}$ . In this case, we can write

(3.24a,b)\begin{equation} \frac{\partial{H_{gy}}}{\partial{p_{\|}}} = \frac{p_{\|}}{m_a \gamma} - \frac{e_a A_{\|}}{m_a c \gamma^3} \left(1 + \frac{2\mu B}{m_a c^2}\right),\quad \boldsymbol{\nabla} H_{gy} = e_a \boldsymbol{\nabla} \langle \psi \rangle,\quad \psi = \phi - \frac{p_{\|}}{m_a c \gamma} A_{\|}. \end{equation}

In a slab, the electromagnetic relativistic gyrokinetic equation takes the form

(3.25)\begin{equation} \frac{\partial{f_{1a}}}{\partial{t}} + \frac{p_{\|}}{m_a \gamma} {\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{1a} = e_a \left({\boldsymbol{b}} \frac{\partial{F_{0a}}}{\partial{p_{\|}}} + c\frac{{\boldsymbol{b}}\times\boldsymbol{\nabla} F_{0a}}{e_a B}\right) \boldsymbol{\cdot}\boldsymbol{\nabla} \langle \psi \rangle. \end{equation}

The relativistic parallel Ampere's law reads

(3.26)\begin{equation} -\frac{1}{4{\rm \pi}} \nabla_{{\perp}}^2 A_{\|} = \sum_a e_a \int {\rm d}^6 Z \delta^3_{gy} \frac{p_{\|}}{m_a c \gamma} \left[f_{1a} + \{S_1, F_{0a}\} + \frac{e_a A_{\|}}{c} {\boldsymbol{b}}\boldsymbol{\cdot}\{{\boldsymbol{R}} + {\boldsymbol{\rho}}, F_{0a}\}\right]. \end{equation}

Evaluating the Poisson brackets (Brizard & Chan Reference Brizard and Chan1999) in the usual way, we obtain

(3.27)\begin{equation} f_{1a} + \{S_1, F_{0a}\} + \frac{e_a A_{\|}}{c} {\boldsymbol{b}} \boldsymbol{\cdot} \{{\boldsymbol{R}} + {\boldsymbol{\rho}}, F_{0a}\} = f_{1a} + \frac{e_a \gamma}{B} \tilde{\psi} \frac{\partial{F_{0a}}}{\partial{\mu}} + \frac{e_a A_{\|}}{c} \frac{\partial{F_{0a}}}{\partial{p_{\|}}}. \end{equation}

Substituting the partial derivatives of the Maxwell–Jüttner distribution function, computed in (3.11a,b) and (3.12), results in

(3.28)\begin{equation} \frac{e_a \gamma}{B} \tilde{\psi} \frac{\partial{F_{0a}}}{\partial{\mu}} + \frac{e_a A_{\|}}{c} \frac{\partial{F_{0a}}}{\partial{p_{\|}}} ={-} \frac{e_a \tilde{\phi}}{T_a} F_{0a} - \frac{p_{\|}}{m_a c \gamma} \frac{e_a \langle A_{\|} \rangle}{T_a} F_{0a}. \end{equation}

Finally, the parallel Ampere's law for relativistic plasmas takes the form

(3.29)\begin{gather} \sum_a \frac{4{\rm \pi} e_a^2}{m_a c^2} \int {\rm d}^6 Z \,\delta^3_{gy} \left(\frac{p_{\|}}{m_a \gamma}\right)^2\frac{m_a F_{0a}}{T_a} \langle A_{\|} \rangle - \nabla_{{\perp}}^2 A_{\|} = \frac{4{\rm \pi}}{c} \sum_a \bar{j}_{1\|a}, \end{gather}

(3.30)\begin{gather}\bar{j}_{1\|a} = e_a \int {\rm d}^6 Z \,\delta^3_{gy} \frac{p_{\|}}{m_a \gamma} f_{1a}. \end{gather}

It is straightforward to combine this equation with the relativistic quasineutrality condition and the gyrokinetic equation in the local limit in order to obtain the relativistic generalization of the electromagnetic dispersion relation (2.13).