Gyrokinetic applications in electron–positron and non-neutral plasmas

Abstract

In this paper, we summarize our recent work on gyrokinetic applications in electron–positron and non-neutral plasma. The electrostatic stability of electron–positron plasmas was investigated in dipole and slab geometries, with and without ion admixture. The gyrokinetic dispersion relation was derived and, for the slab case, extended to non-neutral plasmas. Here, we further extend the gyrokinetic formulation to the relativistic regime. Electron–positron plasmas are found to be remarkably stable as long as perfect symmetry between the two species prevails, but instabilities appear if this symmetry is broken, for instance by the introduction of impurities or magnetic curvature.

1 Introduction

Natural electron–positron plasmas can be found in many places in the Universe. Normally, these are highly energetic locations, such as pulsar magneto-spheres (Spitkovsky 2008) or Poynting-flux dominated astrophysical jets (Lyutikov & Blackman 2001), since photons with very high energies or very strong electromagnetic fields are needed for the pair creation. Compact objects, such as pulsars, can provide the energy sufficient for the pair creation, due to their large masses, extremely fast rotations and very strong magnetic fields. Charged particles drifting in such strong inhomogeneous magnetic fields can provide high-energy gamma radiation via bremsstrahlung. Astrophysical pair plasmas are normally not ‘pure’ and contain other species (Pétri 2016), e.g. protons or iron ions. These plasmas are normally relativistic and may coexist with strong radiation (Uzdensky 2016; Cruz et al. 2021). Charge neutrality of the plasmas surrounding magnetized rotating compact objects (e.g. pulsars) is often violated (Pétri 2009).

Recently, there has been great interest in the production of pair plasmas in a laboratory. Relativistic pair plasmas with properties similar to the astrophysical ones can be obtained in laser experiments where pairs are produced via an interaction of laser beams with a gold target (Chen & Fiuza 2023). In contrast, magnetically confined pair plasmas (Stoneking et al. 2020) will employ external sources of positrons and therefore have low temperatures (few eV). These plasmas are non-relativistic and standard magnetic fusion plasma theory and numerical tools can readily be applied to them. They can be confined in magnetic dipole traps or stellarators and are not necessarily quasineutral (since plasmas with an arbitrary degree of charge neutrality can be confined in an external magnetic field). A well-established application of laboratory non-neutral plasmas with a large number of positrons (Surko & Greaves 2004) is the field of antimatter research which includes positronium (Cassidy & Mills 2007) and antihydrogen (Fajans & Surko 2020) production. An important technique in the antihydrogen creation is the cooling of antiprotons via interactions with an electron gas (Rolston & Gabrielse 1989). Stability of such complex many-component plasmas may also be addressed borrowing well-established tools from the magnetic-confinement research.

Gyrokinetic theory (Brizard & Hahm 2007; Catto 2019) is a reduced description of the low-frequency dynamics in magnetized plasmas. It is a standard model used for turbulence (Garbet et al. 2010) and energetic-particle-driven instabilities (Chen & Zonca 2016) in magnetic fusion research. Starting around 2009, gyrokinetic theory has also been applied to astrophysical systems (Schekochihin et al. 2009), such as the solar wind and accretion disks. A considerable amount of analytical work and dozens of numerical codes exist that employ gyrokinetic theory in magnetic fusion and astrophysical-plasma contexts. It seems promising to apply gyrokinetic theory and numerical tools also for antimatter (electron–positron) and non-neutral (e.g. electron–antiproton) plasma problems. An extension of gyrokinetic theory to the relativistic regimes would be desirable for the astrophysical applications involving compact objects and for the laboratory laser plasmas.

It has been shown (Stenson et al. 2017) that the wave dynamics drastically simplifies in pair plasmas in the cold-plasma limit. A natural question arises as to whether this strong simplification still holds in the gyrokinetic regime. Since the confinement theorem (Dubin & O'Neil 1999) does not apply to quasineutral plasmas, a cylindrical configuration, such as the Penning trap, is not an option for confinement of plasma containing both positively and negatively charged particles. Toroidal configurations are required such as dipoles (Saitoh et al. 2014) or stellarators (Pedersen et al. 2012). Both these options are being currently pursued in ongoing laboratory projects (Saitoh et al. 2014; Stenson 2019). The toroidal geometry violates pair-plasma symmetry since the curvature drift direction depends on the sign of the particle charge. This may lead to collective instabilities, as has recently been shown for the gyrokinetic regime (Mishchenko, Plunk & Helander 2018a). Naturally, these collective micro-instabilities can lead to particle and energy turbulent transport which can be harmful or beneficial depending on the problem at hand. The role of the magnetic-field configuration in pair-plasma stability and resulting confinement is an issue of practical relevance for the experiments under construction (Saitoh et al. 2014; Stenson 2019). Our recent work on gyrokinetic applications in electron–positron and non-neutral plasmas includes the following:

1. (i) The investigation of the electrostatic stability of electron–positron plasmas (Mishchenko et al. 2018a) in a dipole geometry. Here, the kinetic dispersion relation for sub-bounce-frequency instabilities has been derived and solved. For the zero-Debye-length case, the stability diagram has been found to exhibit singular behaviour. However, when the Debye length is non-zero, a fluid mode appears, resolving the observed singularity. It has been demonstrated that both the temperature and density gradients can drive instability.

2. (ii) The study of the gyrokinetic stability of electron–positron plasmas contaminated by ion (proton) admixture (Mishchenko et al. 2018b) in a slab geometry. The appropriate dispersion relation has been derived and solved. The ion-temperature-gradient-driven instability (ITG), the electron-temperature-gradient-driven instability (ETG), the universal mode and the shear Alfvén wave were considered.

3. (iii) The investigation of the confining properties of dipole and stellarator geometries, ranging from pure electron plasmas through to quasineutral. We have shown (Kennedy & Mishchenko 2019) that non-neutral plasmas can be unstable with respect to both density-gradient- and temperature-gradient-driven instabilities.

4. (iv) The numerical study of the gyrokinetic stability of plasmas in different magnetic geometries (Kennedy et al. 2020). The stability of plasmas has been examined varying the mass ratio between the positive and negative charge carriers, from conventional hydrogen plasmas through to electron–positron plasmas. Stability was studied for prescribed temperature and density gradients in an axisymmetric tokamak and a non-axisymmetric quasi-isodynamic stellarator configurations.

5. (v) The linear gyrokinetic simulations of magnetically confined electron–positron plasmas have been performed (Kennedy et al. 2018) in the dipole geometry and parameter regimes likely to be relevant for upcoming laboratory experiments (Stoneking et al. 2020). Our results have demonstrated the existence of unstable entropy modes and interchange modes in pair plasmas.

The paper is organized as follows. In § 2, we review the non-relativistic gyrokinetic theory for pair and non-neutral plasmas. In § 3, the relativistic extension is derived. In § 4, we draw our conclusions.

2 Non-relativistic case

2.1 Slab geometry

Following Helander (2014) and Helander & Connor (2016), 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. (2018b). 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. 2012). These results extend the finding of Helander (2014) 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 2019) 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. (2018a) 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. (2018a). In contrast, one needs some violation of the symmetry between the species for instability in a slab.

3 Relativistic plasmas

3.1 Relativistic gyrokinetic equation

The relativistic gyrokinetic equation has been derived in Brizard & Chan (1999) 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 (1999), 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 (1911), Zenitani (2015) and Cercignani & Kremer (2012),

(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 (1999), 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 1987; Castejón & Pavlov 2006) 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 1999) 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).

4 Conclusions

In this paper, we have summarized our recent work on gyrokinetic applications in electron–positron and non-neutral plasma (Mishchenko et al. 2018a; Kennedy & Mishchenko 2019). The gyrokinetic stability of electron–positron plasmas contaminated by ion (proton) admixture has been studied in Mishchenko et al. (2018b) in a slab geometry. The appropriate dispersion relation was derived and solved. The destabilization of ITGs, ETGs and universal modes at finite ion contamination were considered. It has been shown in Kennedy & Mishchenko (2019) that drift instabilities can be excited in non-neutral plasmas. In Mishchenko et al. (2018a), the electrostatic stability of electron–positron plasmas has been investigated in the dipole geometry. Linear gyrokinetic simulations of magnetically confined electron–positron plasmas were performed in dipole (Kennedy et al. 2020) and stellarator (Kennedy et al. 2018) geometries.

Similarly to the cold-plasma case (Stenson et al. 2017), a drastic reduction of the unstable solutions is found for pure pair plasma also in the gyrokinetic regime. Pure pair plasmas are gyrokinetically stable for all gradients if the symmetry between the species is not violated. However, any species asymmetry can drive gyrokinetic instabilities. The asymmetry can be caused by a contamination (e.g. with protons or other ion species), plasma non-neutrality or different electron and positron temperature profiles. Finally, magnetic curvature can cause gyrokinetic instabilities of pair plasmas too, since the magnetic drift depends on the sign of the particle charge which is opposite for the electrons and positrons. Hence, toroidally confined pair plasmas (in a dipole or a stellarator trap) can be turbulent. This turbulence may lead to the self-organization process such as the inward pinch or zonal flows.

In this paper, the dispersion relation for the slab case is generalized to the relativistic regime. It is found that the remarkable pair-plasma stability can be extended to relativistic applications with the same limitations of a perfect species symmetry and absence of curvature drifts. Formally, relativistic modifications of the dispersion relation and other basic equations (the gyrokinetic and Maxwell equations) are rather moderate. However, the parallel resonance structure becomes more complex due to the presence of the relativistic Lorentz factor which depends both on the parallel momentum $p_{\|}$ and the magnetic moment $\mu$. Also, the relativistic ambient distribution function (such as the Maxwell–Jüttner distribution function) leads to complications since the usual plasma-dispersion-function formalism cannot be applied. Some analytical progress can be made in the weak-relativistic limit but the calculations become quickly very cumbersome.

One limitation of this extension is that, as in Brizard & Chan (1999), only the special relativity contributions have been included in this paper. However, for astrophysical applications in areas of extreme gravitation, such as the surroundings of black holes, general relativity must be used. The appropriate extension of the gyrokinetic theory is being elaborated (Beklemishev & Tessarotto 2004). Another limitation is in the Maxwell–Jüttner distribution function employed in the paper. It assumes some agent, such as Coulomb collisions, establishing a local thermodynamic equilibrium and neglects pair creation and annihilation, which are processes inherent to electron–positron plasmas at high enough energy density. In future, a formulation based on the Dirac equation has to be developed in order to account for these complications (Uzdensky & Rightley 2014). Finally, radiation and its reaction on plasma particles may play an important role, for example affecting the reconnection dynamics (Uzdensky 2016). Such effects seem to be out of scope of the regular gyrokinetic theory because of its low-frequency ordering.