2.1 Gyrokinetics

As is often the case in the study of magnetically confined plasmas, we shall assume that the fluctuations within our plasma obey the standard gyrokinetic ordering (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013 or Catto Reference Catto2019); that is, for fluctuations with a characteristic frequency $\omega$ and wavenumbers $k_\parallel$ and $k_\perp$ parallel and perpendicular to the equilibrium magnetic field direction $\boldsymbol {b}_0 = \boldsymbol {B}_0/B_0$, we have

(2.1)\begin{equation} \frac{\omega}{\varOmega_s} \sim \frac{\nu_{s s'}}{\varOmega_s} \sim \frac{k_\parallel}{k_\perp} \sim \frac{q_s \phi}{T_{0s}} \sim \frac{\delta \! B_\parallel}{B_0} \sim \frac{\delta \! \boldsymbol{B}_\perp}{B_0} \sim \frac{\rho_s}{L} \ll 1, \end{equation}

where $\varOmega _s = q_s B_0/m_s c$ is the cyclotron frequency of species $s$ with charge $q_s$, equilibrium density and temperature $n_{0s}$ and $T_{0s}$, respectively, mass $m_s$ and thermal speed $v_{{\rm th}s} = \sqrt {2T_{0s}/m_s}$, $\nu _{s s'}$ is the typical collision frequency, $\rho _s = v_{{\rm th}s}/|\varOmega _s|$ is the thermal Larmor radius, $\delta \! B_\parallel$ and $\delta \! \boldsymbol {B}_\perp$ are the fluctuations of the magnetic field parallel and perpendicular to the equilibrium direction, respectively, and $L$ is a typical equilibrium length scale. It is assumed that all equilibrium quantities evolve on the (long) transport timescale $\tau _E^{-1} \sim (\rho _s/L)^3 \varOmega _s$, and so will be considered static throughout the remainder of this paper.

Under the ordering of (2.1), the perturbed distribution function $\delta \! f_s$ consists of the Boltzmann and gyrokinetic parts:

(2.2)\begin{equation} \delta \! f_s(\boldsymbol{r},\boldsymbol{v},t) ={-} \frac{q_s \phi(\boldsymbol{r},t)}{T_{0s}} f_{0s}(x,\boldsymbol{v}) + h_s(\boldsymbol{R}_s, v_\perp, v_\parallel,t), \end{equation}

where $\boldsymbol {R}_s = \boldsymbol {r} - \boldsymbol {b}_0 \times \boldsymbol {v}_\perp /\varOmega _s$ is the guiding-centre position, and $h_s$ evolves according to the gyrokinetic equation

(2.3)\begin{equation} \frac{\partial}{\partial t} \left( h_s - \frac{q_s \left\langle \chi \right\rangle_{\boldsymbol{R}_s}}{T_{0s}} f_{0s} \right) + \left(v_\parallel \boldsymbol{b}_0 + \boldsymbol{v}_{d s} \right) \boldsymbol{\cdot} {\boldsymbol{\nabla}} h_s + \boldsymbol{v}_\chi \boldsymbol{\cdot} {\boldsymbol{\nabla}}_\perp \left( h_s + f_{0s} \right) = \left( \frac{\partial h_s}{\partial t} \right)_c. \end{equation}

In the above, and throughout this paper, $\left \langle \cdots \right \rangle _{\boldsymbol {R}_s}$ denotes the standard gyroaverage at constant $\boldsymbol {R}_s$. Here, $\chi = \phi - \boldsymbol {v}\boldsymbol {\cdot } \boldsymbol {A}/c$ is the gyrokinetic potential ($\phi$ and $\boldsymbol {A}$ are the scalar and vector potential, respectively, under the Coulomb gauge ${\boldsymbol {\nabla }} \boldsymbol {\cdot } \boldsymbol {A} = 0$) that gives rise to the drift velocity

(2.4)\begin{equation} \boldsymbol{v}_\chi =\frac{c}{B_0}\boldsymbol{b}_0 \times \frac{\partial \left\langle \chi \right\rangle_{\boldsymbol{R}_s}}{\partial \boldsymbol{R}_s}, \end{equation}

which includes the $\boldsymbol {E} \times \boldsymbol {B}$ drift, the parallel streaming along perturbed field lines and the ${\boldsymbol {\nabla }} \! B$ drift associated with the perturbed magnetic field. This gives rise to nonlinearities (with which we will not be concerned in this paper), as well as the familiar gyrokinetic drive associated with the equilibrium distribution $f_{0s}$, viz.

(2.5)\begin{equation} \boldsymbol{v}_\chi \boldsymbol{\cdot} {\boldsymbol{\nabla}}_\perp f_{0s} ={-} \frac{c}{B_0} \left( \boldsymbol{b}_0 \times \frac{\partial \left\langle \chi \right\rangle_{\boldsymbol{R}_s}}{\partial \boldsymbol{R}_s} \right) \boldsymbol{\cdot} {\boldsymbol{\nabla}} x \left[\frac{1}{L_{n_s}} + \frac{\eta_s}{L_{n_s}} \left( \frac{v^2}{v_{{\rm th}s}^2} - \frac{3}{2} \right) \right] f_{0s}, \end{equation}

where

(2.6)\begin{equation} L_{n_s}^{{-}1} ={-} \frac{1}{n_{0s}} \frac{\partial n_{0s}}{\partial x}, \quad L_{T_s}^{{-}1} ={-} \frac{1}{T_{0s}} \frac{\partial T_{0s}}{\partial x}, \quad \eta_s = \frac{L_{n_s}}{L_{T_s}} \end{equation}

are the characteristic length scales associated with the radial equilibrium gradients of both density and temperature, respectively, $\eta _s$ is their ratio, and $x$ is the direction of the equilibrium gradients. The magnetic drifts associated with the equilibrium field are

(2.7)\begin{equation} \boldsymbol{v}_{ds} = \frac{\boldsymbol{b}_0}{\varOmega_s} \times \left[ v_\parallel^2 \boldsymbol{b}_0\boldsymbol{\cdot} {\boldsymbol{\nabla}}\boldsymbol{b}_0 + \frac{1}{2}v_\perp^2 {\boldsymbol{\nabla}}\log B_0 \right]. \end{equation}

The last term on the right-hand side of (2.3) is the (linearised) collision operator, which we henceforth neglect given that we are interested in studying collisionless dynamics. The electromagnetic fields appearing in the gyrokinetic equation (2.3) are determined by the quasineutrality condition

(2.8)\begin{equation} 0 = \sum_{s} q_s \delta n_s = \sum_s q_s \left[ -\frac{q_s \phi}{T_{0s}} n_{0s} + \int \,\mathrm{d}^3 \boldsymbol{v} \left\langle h_s \right\rangle_{\boldsymbol{r}}\right], \end{equation}

where $\left \langle \cdots \right \rangle _{\boldsymbol {r}}$ denotes the gyroaverage at constant $\boldsymbol {r}$, and by the parallel and perpendicular parts of Ampère's law, which are respectively

(2.9)\begin{align} {\boldsymbol{\nabla}}_\perp^2 A_\parallel&={-} \frac{4{\rm \pi}}{c} \sum_s q_s \int \,\mathrm{d}^3 \boldsymbol{v} \: v_\parallel \left\langle h_s \right\rangle_{\boldsymbol{r}}, \end{align}

(2.10)\begin{align}{\boldsymbol{\nabla}}_\perp^2 \delta \! B_\parallel&={-} \frac{4{\rm \pi}}{c} \boldsymbol{b}_0 \boldsymbol{\cdot} \left[ {\boldsymbol{\nabla}}_\perp{\times} \sum_s q_s \int \,\mathrm{d}^3 \boldsymbol{v} \left\langle \boldsymbol{v}_\perp h_s \right\rangle_{\boldsymbol{r}} \right]. \end{align}

Together, (2.3) and (2.8)–(2.10) form a closed system of equations that, in principle, allows us to determine $h_s$ and thus the evolution of the fluctuations in our plasma. In this work, we solve the linear part of this system in the ‘local’ limit (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995): we assume that the gradients of all equilibrium quantities are constant – including the geometric coefficients $\boldsymbol {b}_0\boldsymbol {\cdot }{\boldsymbol {\nabla }}\boldsymbol {b}_0$ and ${\boldsymbol {\nabla }}\log B_0$ that appear in the magnetic drifts – and choose orthonormal coordinates $(x, y, z)$, in which $\hat {\boldsymbol {z}}=\boldsymbol {b}_0$ is the direction of the magnetic field, $\hat {\boldsymbol {x}}$ is, as above, the direction of the equilibrium gradients (cf. the radial direction in toroidal geometry), and $\hat {\boldsymbol {y}} \equiv \boldsymbol {b}_0 \times \hat {\boldsymbol {x}}$ is the binormal direction (cf. the poloidal direction in toroidal geometry). One can think of this geometry as that of a $Z$-pinch (see Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022; Ivanov, Schekochihin & Dorland Reference Ivanov, Schekochihin and Dorland2022) due to the assumption of constant magnetic curvature and lack of magnetic shear, which we have implicitly assumed. Under these assumptions, the system of (2.3), (2.8)–(2.10) is homogeneous in space, allowing us to impose periodic boundary conditions in all three spatial dimensions.

In the next section, we consider the time evolution of a single Fourier mode and obtain the resulting gyrokinetic dispersion relation.

2.2 Linear gyrokinetic problem

Neglecting the nonlinear term and introducing the spatial Fourier decomposition:

(2.11)\begin{equation} h_s(\boldsymbol{R}_s, v_\perp, v_\parallel, t) = \sum_{\boldsymbol{k}} {h_s}_{\boldsymbol{k}}(v_\perp, v_\parallel, t) {\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{R}_{s}}, \quad \chi(\boldsymbol{r},t) = \sum_{\boldsymbol{k}}{\chi}_{\boldsymbol{k}}(t) {\rm e}^{{\rm i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r} }, \end{equation}

with $\boldsymbol {k} = \boldsymbol {k_\perp } + k_\parallel \boldsymbol {b}_0$, the Fourier modes ${h_s}_{\boldsymbol {k}}$ and ${\chi }_{\boldsymbol {k}}$ can be shown to satisfy

(2.12)\begin{equation} \frac{\partial}{\partial t} \left( {h_s}_{\boldsymbol{k}} - \frac{q_{s} \left\langle{\chi}_{\boldsymbol{k}} \right\rangle_{\boldsymbol{R}_{s}}}{T_{0s}}f_{0s} \right) + {\rm i} k_\parallel v_\parallel {h_s}_{\boldsymbol{k}} + {\rm i} \omega_{Ds} {h_s}_{\boldsymbol{k}} - {\rm i} \omega_{*s}^T \frac{q_{s} \left\langle{\chi}_{\boldsymbol{k}} \right\rangle_{\boldsymbol{R}_{s}}}{T_{0{s}}}f_{0{s}} = 0, \end{equation}

where we have defined the drift frequencies associated with the equilibrium gradients of species $s$ [cf. (2.5)]:

(2.13)\begin{equation} \omega_{*s}^T = \omega_{*s} \left[ 1 + \eta_s \left( \frac{v^2}{v_{{\rm th}s}^2} - \frac{3}{2} \right) \right], \quad \omega_{*s} ={-} \frac{k_y c T_{0s}}{q_s B_0 L_{n_s}}, \end{equation}

and with the equilibrium magnetic field curvature and gradient [cf. (2.7)]:

(2.14)\begin{equation} \omega_{Ds} = \frac{2v_\parallel^2}{v_{{\rm th}s}^2} \omega_{\kappa s} + \frac{v_\perp^2}{v_{{\rm th}s}^2} \omega_{{\boldsymbol{\nabla}} \! B s}, \end{equation}

where

(2.15)\begin{equation} \omega_{\kappa s} = \frac{v_{{\rm th}s}^2 }{2\varOmega_s} \boldsymbol{k}_\perp \boldsymbol{\cdot} \left[\boldsymbol{b}_0 \times (\boldsymbol{b}_0 \boldsymbol{\cdot} {\boldsymbol{\nabla}} )\boldsymbol{b}_0 \right], \quad \omega_{{\boldsymbol{\nabla}} \! Bs} = \frac{v_{{\rm th}s}^2}{2\varOmega_s} \boldsymbol{k}_\perp \boldsymbol{\cdot} \left(\boldsymbol{b}_0 \times {\boldsymbol{\nabla}} \log B_0 \right). \end{equation}

Starting from the perpendicular force balance of the gyrokinetic equilibrium [see (128) in Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013], it is straightforward to show that the difference between these two drifts is given by

(2.16)\begin{equation} \omega_{\kappa s} - \omega_{{\boldsymbol{\nabla}} \! B s} = \frac{v_{{\rm th}s}^2}{2\varOmega_s} \boldsymbol{k_\perp} \boldsymbol{\cdot} \left( \boldsymbol{b}_0 \times {\boldsymbol{\nabla}} x \right) \left.\frac{\partial}{\partial x}\right|_{B_0} \sum_{s'} \frac{\beta_{s'}}{2} , \end{equation}

where $\beta _s = 8 {\rm \pi}n_{0s} T_{0s}/B_0^2$ is the plasma beta of species $s$. Lastly, the gyroaveraged Fourier-transformed gyrokinetic potential is

(2.17)\begin{equation} \left\langle {\chi}_{\boldsymbol{k}} \right\rangle_{\boldsymbol{R}_s} = {\rm J}_0(b_s) \left({\phi}_{\boldsymbol{k}} - \frac{v_\parallel A_{{\parallel} \boldsymbol{k}}}{c} \right) + \frac{2 {\rm J}_1 (b_s)}{b_s} \frac{T_{0s}}{q_s} \frac{v_\perp^2 }{v_{{\rm th}s}^2}\frac{{\delta \! B}_{{\parallel} \boldsymbol{k}}}{B_0}, \end{equation}

while the field equations (2.8)–(2.10) can be written as (see, e.g. Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006)

(2.18)\begin{align} \sum_s \frac{q_s^2 n_{0s}}{T_{0s}} {\phi}_{\boldsymbol{k}} &= \sum_s q_s \int \,\mathrm{d}^3 \boldsymbol{v} \: {\rm J}_0(b_s) \tilde{h}_s, \end{align}

(2.19)\begin{align}k_\perp^2 A_{{\parallel} \boldsymbol{k}} &= \frac{4{\rm \pi}}{c} \sum_s q_s \int \,\mathrm{d}^3 \boldsymbol{v} \: v_\parallel {\rm J}_0(b_s)\tilde{h}_s , \end{align}

(2.20)\begin{align}\frac{{\delta \! B}_{{\parallel} \boldsymbol{k}}}{B_0} &={-} \frac{1}{2} \sum_s \frac{\beta_s}{n_{0s}} \int \,\mathrm{d}^3 \boldsymbol{v} \: \frac{v_\perp^2}{v_{{\rm th}s}^2} \frac{2 {\rm J}_1 (b_s)}{b_s} \tilde{h}_s, \end{align}

where $b_s = k_\perp v_\perp /\varOmega _s$, and ${\rm J}_0$, ${\rm J}_1$ are the Bessel functions of the first kind (Abramowitz & Stegun Reference Abramowitz and Stegun1972) that capture finite-Larmor-radius effects. It will prove convenient to combine ${\phi }_{\boldsymbol {k}}$, $A_{\parallel \boldsymbol {k}}$, and ${\delta \! B}_{\parallel \boldsymbol {k}}$ into a single vector $\boldsymbol {\chi }_{\boldsymbol {k}}$ given by

(2.21)\begin{equation} \boldsymbol{\chi}_{\boldsymbol{k}} = \left( \frac{q_r {\phi}_{\boldsymbol{k}}}{T_{0r}}, \: \frac{k_\parallel }{|k_\parallel| } \frac{A_{{\parallel} \boldsymbol{k}}}{\rho_r B_0}, \: \frac{{\delta \! B}_{{\parallel} \boldsymbol{k}}}{B_0} \right)^{\rm T}. \end{equation}

Here, and in what follows, we normalise the electromagnetic fields using an arbitrary reference mass $m_r$, density $n_{0r}$, thermal velocity $v_{{\rm th}r}$, temperature $T_{0r}$, and gyroradius $\rho _r$.

Following Landau (Reference Landau1946), we consider an initial-value problem and introduce the Laplace transformations

(2.22)\begin{equation} {\hat{h}}_{s \boldsymbol{k}}(v_\perp, v_\parallel, p) = \int_0^{\infty} \,\mathrm{d} t \; {\rm e}^{{-}pt} {h_s}_{\boldsymbol{k}}(v_\perp, v_\parallel, t), \quad \hat{\boldsymbol{\chi}}_{\boldsymbol{k}}(p) = \int_0^{\infty} \,\mathrm{d} t \; {\rm e}^{{-}pt} \boldsymbol{\chi}_{\boldsymbol{k}}(t). \end{equation}

Assuming there exist positive real $m$ and $M$ such that

(2.23)\begin{equation} \left|{h_s}_{\boldsymbol{k}}(v_\perp, v_\parallel, t)\right|, \quad \left| \boldsymbol{\chi}_{\boldsymbol{k}}(t) \right| \leqslant Me^{m t}, \end{equation}

for all $t > 0$, and picking any real $\sigma$ with $\sigma > m$, the integrals in (2.22) converge and the transformed distributions ${\hat {h}}_{s \boldsymbol {k}}$ and fields $\hat {\boldsymbol {\chi }}_{\boldsymbol {k}}$ are analytic for all complex values of $p$ with $\mathrm {Re}(p) \geqslant \sigma$. The inverse transformations are given by

(2.24)\begin{gather} {h_s}_{\boldsymbol{k}}(v_\perp, v_\parallel, t) = \frac{1}{2 {\rm \pi}{\rm i}} \int_{{C_\sigma}} \,\mathrm{d} p \: {\rm e}^{pt} {\hat{h}}_{s \boldsymbol{k}}(v_\perp, v_\parallel, p), \quad \boldsymbol{\chi}_{\boldsymbol{k}}(t) = \frac{1}{2 {\rm \pi}{\rm i}} \int_{{C_\sigma}} \,\mathrm{d} p \: {\rm e}^{pt} \hat{\boldsymbol{\chi}}_{\boldsymbol{k}}(p), \end{gather}

where the contour of integration ${C_\sigma }$ is along a straight line parallel to the imaginary axis and intersecting the real axis at $\mathrm {Re}(p) = \sigma$, as in figure 1 (this is the so-called Bromwich contour).

Figure 1.

The complex $p$ plane, with $\mathrm {Re}(p)$ and $\mathrm {Im}(p)$ shown on the horizontal and vertical axes, respectively. The contour of integration for the inverse Laplace transform ${C_\sigma }$ is is a vertical straight line at $\mathrm {Re}(p) = \sigma$, to the right of which (i.e. in the shaded grey region) the functions ${\hat {h}}_{s \boldsymbol {k}}$ and ${\hat {\chi }}_{\boldsymbol {k}}$ are guaranteed to be analytic. Singularities, such as poles (indicated by crosses) or branch cuts (indicated by the zigzag line), could exist at $\mathrm {Re}(p) < \sigma$.

Performing the Laplace transform as in (2.22), (2.12) straightforwardly becomes

(2.25)\begin{equation} {\hat{h}}_{s \boldsymbol{k}} = \frac{p + {\rm i} \omega_{*s}^T}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}} \frac{q_{s} \left\langle{\hat{\chi}}_{\boldsymbol{k}} \right\rangle_{\boldsymbol{R}_{s}}}{T_{0s}}f_{0s} + \frac{{{g}_{s \boldsymbol{k}}}}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}}, \end{equation}

where ${g}_{s \boldsymbol {k}}$ is the initial condition:

(2.26)\begin{equation} {g}_{s \boldsymbol{k}}(v_\perp, v_\parallel) = {h_s}_{\boldsymbol{k}}(v_\perp, v_\parallel, {t=0}) - \frac{q_{s} \left\langle{\hat{\chi}}_{\boldsymbol{k}} ({t=0}) \right\rangle_{\boldsymbol{R}_{s}}}{T_{0s}}f_{0s}. \end{equation}

Then, normalising the characteristic frequencies to the parallel-streaming rateFootnote ¹

(2.27)\begin{equation} \zeta_s = \frac{ip}{|k_\parallel| v_{{\rm th}s}}, \quad \zeta_{*s} = \frac{\omega_{*s}}{|k_\parallel| v_{{\rm th}s}}, \quad \zeta_{\kappa s} = \frac{\omega_{\kappa s}}{|k_\parallel| v_{{\rm th}s}}, \quad \zeta_{{\boldsymbol{\nabla}} \! B s} = \frac{\omega_{{\boldsymbol{\nabla}} \! B s}}{|k_\parallel| v_{{\rm th}s}}, \end{equation}

and defining the dimensionless velocity variables

(2.28)\begin{equation} u = \frac{k_\parallel }{|k_\parallel| } \frac{v_\parallel}{v_{{\rm th}s}}, \quad \mu = \frac{v_\perp^2}{v_{{\rm th}s}^2}, \end{equation}

we substitute (2.25) into the Laplace transforms of the field equations (2.18)–(2.20) to obtain the linear eigenvalue problem

(2.29)\begin{equation} \boldsymbol{\mathsf{L}} \hat{\boldsymbol{\chi}}_{\boldsymbol{k}} + \boldsymbol{G} = 0 , \end{equation}

in which $\boldsymbol{\mathsf{L}} $ is the linear coefficient matrix and $\boldsymbol {G}$ is the vector of the initial conditions of the fields. The components of $\boldsymbol{\mathsf{L}} $ are given by

(2.30)\begin{align} {\mathsf{L}}_{\phi \phi} &={-} \sum_s \frac{q_s^2 n_{0s} T_{0r}}{q_r^2 n_{0r} T_{0s}} \left\{ 1+ \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \left.\mathcal{I}_{a,b}^{(s)} \right|_{a=b=1}\right\} , \end{align}

(2.31)\begin{align}{\mathsf{L}}_{\phi A} &= 2\sum_s \frac{q_s^2 n_{0s} v_{{\rm th}s} T_{0r}}{q_r^2 n_{0r} v_{{\rm th}r} T_{0s}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \left.\mathcal{J}_{a,b}^{(s)}\right|_{a=b=1}, \end{align}

(2.32)\begin{align}{\mathsf{L}}_{\phi B} &= \sum_s \frac{q_s n_{0s}}{q_r n_{0r}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_b \left.\mathcal{K}_{a,b}^{(s)}\right|_{a=b=1}, \end{align}

(2.33)\begin{align}{\mathsf{L}}_{A \phi} &={-} \sum_s \frac{q_s^2 n_{0s} v_{{\rm th}s} T_{0r}}{q_r^2 n_{0r} v_{{\rm th}r} T_{0s}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \left.\mathcal{J}_{a,b}^{(s)}\right|_{a=b=1} , \end{align}

(2.34)\begin{align}{\mathsf{L}}_{AA} &={-} \frac{B_0^2 ( k_\perp \rho_r)^2}{8 {\rm \pi}n_{0r} T_{0r}} - 2 \sum_s \frac{q_s^2 n_{0s} m_r}{q_r^2 n_{0r} m_s} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_a \left.\mathcal{I}_{a,b}^{(s)}\right|_{a=b=1}, \end{align}

(2.35)\begin{align}{\mathsf{L}}_{AB} &= \sum_s \frac{q_s n_{0s} v_{{\rm th}s}}{q_r n_{0r} v_{{\rm th}r}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_b \left.\mathcal{L}_{a,b}^{(s)}\right|_{a=b=1}, \end{align}

(2.36)\begin{align}{\mathsf{L}}_{B\phi} &={-} \sum_s \frac{\beta_s}{2} \frac{q_s T_{0r}}{q_r T_{0s}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_b \left.\mathcal{K}_{a,b}^{(s)}\right|_{a=b=1}, \end{align}

(2.37)\begin{align}{\mathsf{L}}_{B A} &= \sum_s \beta_s \frac{q_s T_{0r} v_{{\rm th}s}}{q_r T_{0s} v_{{\rm th}r}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_b \left.\mathcal{L}_{a,b}^{(s)} \right|_{a=b=1}, \end{align}

(2.38)\begin{align}{\mathsf{L}}_{BB} &={-} 1 + \sum_s \frac{\beta_s}{2} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \partial_b + \frac{3}{2} \right) \right] \partial_b^2 \left.\mathcal{M}_{a,b}^{(s)}\right|_{a=b=1} , \end{align}

where we have defined the following integrals

(2.39)\begin{align} \mathcal{I}_{a,b}^{(s)} &= \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{ {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} {\rm J}_0^2(b_s), \end{align}

(2.40)\begin{align}\mathcal{J}_{a,b}^{(s)} &= \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{ u {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} {\rm J}_0^2(b_s), \end{align}

(2.41)\begin{align}\mathcal{K}_{a,b}^{(s)} &= \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{ {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} \frac{2 {\rm J}_0 (b_s) {\rm J}_1 (b_s)}{b_s}, \end{align}

(2.42)\begin{align}\mathcal{L}_{a,b}^{(s)} &= \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{ u {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} \frac{2 {\rm J}_0 (b_s) {\rm J}_1 (b_s)}{b_s}, \end{align}

(2.43)\begin{align}\mathcal{M}_{a,b}^{(s)} &= \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{ {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} \left[\frac{2 {\rm J}_1 (b_s)}{b_s} \right]^2. \end{align}

Here, and throughout the remainder of this paper, the parameters $a$ and $b$ are assumed to be both real and positive, ensuring integral convergence. Finally, the components of $\boldsymbol {G} = (G_\phi, G_A, G_B)^{\rm T}$ are given by

(2.44)\begin{align} G_\phi &= \sum_s \frac{q_s n_{0s}}{q_r n_{0r}} \frac{1}{n_{0s}} \int \,\mathrm{d}^3 \boldsymbol{v} \: \frac{{g}_{s \boldsymbol{k}}}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}}{\rm J}_0 (b_s), \end{align}

(2.45)\begin{align}G_A &= \sum_s \frac{q_s n_{0s} v_{{\rm th}s}}{q_r n_{0r} v_{{\rm th}r} } \frac{1}{n_{0s}}\int \,\mathrm{d}^3 \boldsymbol{v} \: \frac{v_\parallel}{v_{{\rm th}s}} \frac{{g}_{s \boldsymbol{k}}}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}}{\rm J}_0 (b_s) , \end{align}

(2.46)\begin{align}G_B &={-} \sum_s \frac{\beta_s}{2}\frac{1}{n_{0s}}\int \,\mathrm{d}^3 \boldsymbol{v} \: \frac{v_\perp^2}{v_{{\rm th}s}^2} \frac{{g}_{s \boldsymbol{k}}}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}}\frac{2 {\rm J}_1(b_s)}{b_s}. \end{align}

The eigenvalue problem (2.29) can be inverted to solve for the fields in the usual way, viz.

(2.47)\begin{equation} \hat{\boldsymbol{\chi}}_{\boldsymbol{k}}(p) = \frac{(\text{adj}\:\boldsymbol{\mathsf{L}}) \boldsymbol{G} }{\det\boldsymbol{\mathsf{L}}}, \end{equation}

where $\text {adj}\:\boldsymbol{\mathsf{L}}$ and $\det \boldsymbol{\mathsf{L}}$ are the adjugate matrix and determinant of the linear matrix $\boldsymbol{\mathsf{L}}$, respectively. The time-dependent fields are then determined by the inverse Laplace transform of (2.47). As discussed above, the integrals in (2.24) are, before analytic continuation, defined for $\mathrm {Re}(p) \geqslant \sigma > 0$. For these values of $p$, $\mathrm {Im}(\zeta _s) > 0$, and so the integrals in (2.39)–(2.43) converge and are analytic functions of $p$. Note that the equation

(2.48)\begin{equation} D(p)\equiv\det\boldsymbol{\mathsf{L}}(p) = 0 \end{equation}

is commonly known as the ‘dispersion relation’, while we shall refer to $D$ itself as the ‘dispersion function’.

2.3 Drift-kinetic limit

To evaluate the integrals (2.39)–(2.43), we specialise to the drift-kinetic limit, in which the perpendicular wavenumbers of the perturbations are assumed small in comparison to the species’ gyroradii, viz.

(2.49)\begin{equation} b_s \sim k_\perp \rho_s \ll 1. \end{equation}

In this limit, the Bessel functions can be expanded as

(2.50)\begin{equation} {\rm J}_0(b_s) = 1 + O(b_s^2), \quad 2{\rm J}_1(b_s)/b_s = 1 + O(b_s^2), \end{equation}

meaning that, to leading order in $b_s$, the contributions of the Bessel functions to the integrals (2.39)–(2.43) are equal to one, and we may write

(2.51)\begin{equation} {\rm I}_{a,b}^{(s)} = \mathcal{I}_{a,b}^{(s)} = \mathcal{K}_{a,b}^{(s)} = \mathcal{M}_{a,b}^{(s)} , \quad {\rm J}_{a,b}^{(s)} = \mathcal{J}_{a,b}^{(s)} = \mathcal{L}_{a,b}^{(s)}, \end{equation}

where ${\rm I}_{a,b}^{(s)}$ and ${\rm J}_{a,b}^{(s)}$ are given by

(2.52)\begin{gather} {\rm I}_{a,b}^{(s)}(\zeta_{s}, \zeta_{\kappa s}, \zeta_{Bs}) = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{{\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)}, \end{gather}

(2.53)\begin{gather}{\rm J}_{a,b}^{(s)}(\zeta_{s}, \zeta_{\kappa s}, \zeta_{Bs}) = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{u {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)}. \end{gather}

Furthermore, we consider the particular case in which the difference between the curvature and ${\boldsymbol {\nabla }} \! B$ drifts, given by the right-hand side of (2.16) is zero and so their associated drift frequencies can be taken to be equal, viz.

(2.54)\begin{equation} \omega_{\kappa s} = \omega_{{\boldsymbol{\nabla}} \! B s} \equiv \omega_{ds} \quad \Rightarrow \quad \zeta_{\kappa s} = \zeta_{{\boldsymbol{\nabla}} \! B s} \equiv \zeta_{ds}. \end{equation}

Note that neither approximation should be interpreted as a consequence of some asymptotic ordering of the parameters describing our gyrokinetic system of equations. Instead, they should be viewed as formal approximations that allow us to obtain a solvable case of a more general one. Their relaxation is discussed in § 7.

With these simplifications, we have reduced our problem to the evaluation of

(2.55)\begin{gather} {\rm I}_{a,b}(\zeta, \zeta_d) = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{{\rm e}^{{-}a u^2 - b \mu}}{u - \zeta + \zeta_d\left(2u^2 + \mu \right)}, \end{gather}

(2.56)\begin{gather}{\rm J}_{a,b}(\zeta, \zeta_d) = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu \: \frac{u {\rm e}^{{-}a u^2 - b \mu}}{u - \zeta +\zeta_d \left(2u^2 + \mu \right)}, \end{gather}

where we have used (2.54) and have dropped the species index for the sake of compactness of notation.

3.1 The plasma dispersion function and Landau's solution

In the absence of magnetic drifts (i.e. when $\zeta _d = 0$), (2.55) and (2.56) can straight- forwardly be written in terms of the well-studied plasma dispersion function (Faddeeva & Terent'ev Reference Faddeeva and Terent'ev1954; Fried & Conte Reference Fried and Conte1961):

(3.1)\begin{equation} Z(\zeta) \equiv \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \: \frac{{\rm e}^{{-}u^2}}{u - \zeta}, \end{equation}

where the integral is defined for $\mathrm {Im}(\zeta )>0$ with the integration contour along the real $u$ axis, as in figure 2(a). In particular, we have that

(3.2)\begin{equation} \left.{\rm I}_{a,b} \right|_{\zeta_d=0} = \frac{1}{b}Z_a(\zeta), \quad \left.{\rm J}_{a,b}\right|_{\zeta_d=0} = \frac{1}{b}\left[\frac{1}{\sqrt{a}} + \zeta Z_a(\zeta)\right], \end{equation}

where we have, for the sake of brevity, introduced the shorthand notation

(3.3)\begin{equation} Z_a(\zeta) \equiv Z(\sqrt{a} \zeta). \end{equation}

The integral in (3.1) can be analytically continued to $\mathrm {Im}(\zeta ) \leqslant 0$ by deforming the contour of integration in such a way as to always keep the pole above it, as shown in figure 2(b, c). This is known as the Landau prescription, and the resultant contour is the well-known Landau contour $C_L$ (Landau Reference Landau1946).

Figure 2.

The Landau prescription for the contour of integration $C_L$ that gives the analytic continuation of (3.1). As the Laplace transform demands $\mathrm {Re}(p) \geqslant \sigma > 0$, the pole $u = \zeta$ is located in the upper-half plane [where $\mathrm {Im}(\zeta ) > 0$, see footnote 1], above the contour of integration, as in panel (a). Therefore, the appropriate analytic continuation for $\mathrm {Re}(p) \leqslant 0$ [i.e. $\mathrm {Im}(\zeta ) \leqslant 0$] demands that the contour must be deformed so as to always remain below the pole, as in panels (b, c). Cauchy's integral theorem ensures that we are free to deform the contour without changing the value of the integral, so long as it does not cross the pole.

The plasma dispersion function (3.1) is ubiquitous in calculations of linear waves and instabilities in systems with a spatially uniform magnetic field; notable examples include the electron-temperature-gradient (see, e.g. Liu Reference Liu1971; Lee et al. Reference Lee, Dong, Guzdar and Liu1987) and ion-temperature-gradient (see, e.g. Rudakov & Sagdeev Reference Rudakov and Sagdeev1961; Coppi et al. Reference Coppi, Furth, Rosenbluth and Sagdeev1966; Sauter, Vaclavik & Skiff Reference Sauter, Vaclavik and Skiff1990; Brunner & Vaclavik Reference Brunner and Vaclavik1998; Smolyakov, Yagi & Kishimoto Reference Smolyakov, Yagi and Kishimoto2002) instabilities, the latter of which we shall consider in § 8. It is also worth noting that the Bessel functions can easily be incorporated into the integrals if $\zeta _d=0$ because the resonant denominators are independent of $\mu$. The resulting expressions involve modified Bessel functions and are well known in the literature (see, e.g. Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006).

3.2 Two-dimensional limit

In the 2-D limit, $k_\parallel \rightarrow 0$ with $\zeta \sim \zeta _d \rightarrow \infty$, it can be shown (via, e.g. a partial-fractions expansion of the integrand) that (2.55) can be expressed exactly in terms of products of the plasma dispersion function (Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989), viz.

(3.4)\begin{equation} {\rm I}_{1,1} ={-} \frac{1}{2 \zeta_d } Z(\sqrt{\varOmega})^2 + {O}{\left(\zeta_d^{{-}2}\right)}, \quad {\rm J}_{1,1} = {O}{\left(\zeta_d^{{-}2}\right)}, \quad \varOmega = \frac{\zeta}{2 \zeta_d}, \end{equation}

with the integral for ${\rm J}_{a,b}$ vanishing to leading order because the integrand in (2.56) is manifestly odd in $u$ in this limit. The analytic continuation for (3.4) is significantly more subtle than in the case of the plasma dispersion function (3.1), owing to the presence of the branch point at $\zeta = 0$; we shall delay discussion of these subtleties until § 6. The solution (3.4) has been used extensively in the investigation of 2-D ITG instabilities (see, e.g. Similon et al. Reference Similon, Sedlak, Stotler, Berk, Horton and Choi1984; Biglari et al. Reference Biglari, Diamond and Rosenbluth1989; Kuroda et al. Reference Kuroda, Sugama, Kanno, Okamoto and Horton1998; Sugama Reference Sugama1999; Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006; Helander et al. Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011; Mishchenko, Plunk & Helander Reference Mishchenko, Plunk and Helander2018; Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018).

3.3 Numerical methods

Owing to their analytical complexity, previous literature has also been devoted to the numerical evaluation of (2.55) and (2.56) (see Beer & Hammett Reference Beer and Hammett1996; Gürcan Reference Gürcan2014; Gültekin & Gürcan Reference Gültekin and Gürcan2018; Gültekin & Gürcan Reference Gültekin and Gürcan2020; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020 and references contained therein). In many cases, this involves expressing these integrals in terms of one-dimensional integrals. For example, writing

(3.5)\begin{equation} \frac{1}{u - \zeta + \zeta_d(2u^2 +\mu)} = {\rm i} \int_{0}^{\mathrm{sgn}[\mathrm{Im} (\zeta)] \infty} \,\mathrm{d} \lambda\: {\rm e}^{{-{\rm i}\lambda[u - \zeta +\zeta_d(2u^2 + \mu)]}} \end{equation}

allows the integration over $u$ and $\mu$ in (2.55) and (2.56) to be done analytically, leaving an integral over $\lambda$ that can be evaluated numerically (cf. Beer & Hammett Reference Beer and Hammett1996; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020). While this method is quite general – in that it also allows the direct inclusion of the Bessel functions in (2.39)–(2.43) – the numerical evaluation of the resultant expressions can often be slow, numerical errors may be difficult to quantify, and subtleties like multivaluedness and branch cuts easy to overlook. This motivates the goal of the present study, viz. to find expressions for these integrals in terms of known functions that can be better understood analytically and more readily computed numerically.

4.1 Multivaluedness

To begin, it shall be useful to consider the integral over $u$ separately, and so we write (2.55) as follows:

(4.1)\begin{equation} {\rm I}_{a,b} = \int_0^\infty \,\mathrm{d} \mu \: {\rm e}^{{-}b\mu} \tilde{{\rm I}}_{a}, \quad \tilde{{\rm I}}_{a} = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \frac{{\rm e}^{{-}au^2}}{u - \zeta + \zeta_d(2u^2 + \mu)}. \end{equation}

Now, for each value of $\mu$, the denominator of $\tilde {{\rm I}}_{a}$ has two zeros at

(4.2)\begin{equation} u=\frac{-1 \pm \sqrt{1 + 8\zeta_d(\zeta - \zeta_d\mu)}}{4\zeta_d} \end{equation}

that produce poles on opposite sides of the integration contour along the real $u$ axis. Unsurprisingly, given the presence of square roots in (4.2), $\tilde {{\rm I}}_{a}$ is a multivalued function. In particular, we shall find that $\tilde {{\rm I}}_{a}$, and thus ${\rm I}_{a,b}$, has two branches, just like the square root. To define these two branches, we need to choose a branch cut, which will allow us to ‘label’ the two zeros in (4.2). Note that this choice cannot (and does not) affect the time evolution of the potentials that results from the inverse Laplace transform of (2.21). It turns out to be analytically convenient to consider the ‘principal’ branch cut for the square-root function, for which $\sqrt {z}$ is discontinuous across $\mathrm {Re}(z) < 0$. We can then define the two branches of the square root, $\sqrt [+]{z}$ and $\sqrt [-]{z}$, where the principal branch satisfies $\sqrt [+]{z} > 0$ for all positive real $z$, and $\mathrm {sgn}[\mathrm {Im} (\sqrt [+]{z})] = \mathrm {sgn}[\mathrm {Im}(z)]$.

At this point, it is non-trivial to define the second branch of ${\rm I}_{a,b}$. The choice of a branch for the square root does not determine the branch of the integral (4.1) but only the labels of the zeros in (4.2) – observe that (4.1) makes no reference to any multivalued functions. Indeed, the function ${\rm I}_{a,b}$ is defined as the integral in (4.1) only for ${\mathrm {Im}(\zeta ) > 0}$; the multivaluedness becomes relevant after one considers the analytic continuation to $\mathrm {Im}(\zeta ) < 0$. To make this explicit, until we perform said continuation, we will make use of the labels $\tilde {{\rm I}}_{a}^+$ and ${\rm I}_{a,b}^+$ to indicate that our expressions only apply to this one branch.

Choosing to work with $\sqrt [+]{}$, the zeros in (4.2) can be written as

(4.3)\begin{equation} u ={\mp} u_\pm, \quad u_\pm \equiv \frac{\sqrt[+]{1 + 8\zeta_d(\zeta - \zeta_d\mu)} \pm 1}{4\zeta_d}. \end{equation}

Using a partial-fraction expansion of the integrand, it follows that

(4.4)\begin{equation} \tilde{{\rm I}}_{a}^+= \frac{1}{2\zeta_d(u_+ + u_-)}\left( \frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^\infty \,\mathrm{d} u \frac{{\rm e}^{{-}au^2}}{u - u_-} - \frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^\infty \,\mathrm{d} u \frac{{\rm e}^{{-}au^2}}{u + u_+ }\right). \end{equation}

Now, given that $\mathrm {Im}(\zeta ) > 0$, the sign of the imaginary part of $\sqrt [+]{1 + 8\zeta _d(\zeta - \zeta _d\mu )}$ is determined by the sign of $\zeta _d$, viz.

(4.5)\begin{equation} \mathrm{sgn}\left[\mathrm{Im}{\sqrt[+]{1 + 8\zeta_d(\zeta - \zeta_d\mu)}}\right] = \mathrm{sgn}(\zeta_d), \end{equation}

and so (4.3) implies that $\mathrm {Im}(u_\pm ) > 0$. Therefore, the first integral in the brackets in (4.4) is manifestly the plasma dispersion function, as the imaginary part of the pole at $u = u_-$ has the correct sign for the definition (3.1), i.e. $\mathrm {Im}(u_-) > 0$. The second integral has a pole at $u = -u_+$ with the opposite sign of its imaginary part, i.e. $\mathrm {Im}(-u_+) < 0$, meaning that it can also be turned into a plasma dispersion function under a straightforward change of variables $u \mapsto -u$ (this effectively flips the pole from being below the real $u$ axis to being above it). Thus, it follows that (4.4) can be written as

(4.6)\begin{equation} \tilde{{\rm I}}_{a}^+= \frac{1}{2\zeta_d} \frac{Z_a(u_+) + Z_a(u_- )}{u_+ + u_-}, \end{equation}

where we have used the shorthand notation (3.3) for $Z_a$.

4.2 Explicit evaluation of ${\rm I}_{a,b}^+$

Using (4.6), our expression for ${\rm I}_{a,b}^+$ thus becomes:

(4.7)\begin{equation} {\rm I}_{a,b}^+= \frac{1}{2\zeta_d} \int_0^\infty \,\mathrm{d} \mu \: {\rm e}^{{-}b \mu} \frac{Z_a(u_+) + Z_a(u_- )}{u_+ + u_-}. \end{equation}

Using (4.3), together with the property $Z'(u) = -2 [1 + uZ(u)]$, it can be deduced that

(4.8)\begin{equation} \mu = \frac{\zeta}{\zeta_d} + \frac{1}{4 \zeta_d^2} - \left(u_+^2 + u_-^2 \right) \end{equation}

and

(4.9)\begin{equation} \frac{{\rm e}^{a u_\pm^2}}{u_+ + u_-} = \frac{1}{\sqrt{a}} \frac{\partial}{\partial \mu } \left[{\rm e}^{a u_\pm^2} Z_a(u_\pm) \right]. \end{equation}

It is then a matter of straightforward algebra to show that (4.7) can be rewritten as

(4.10)\begin{equation} {\rm I}_{a,b}^+= \frac{1}{2 \sqrt{a} \zeta_d } \int_0^\infty \,\mathrm{d} \mu \: {\rm e}^{-(b-a) \mu} \frac{\partial}{\partial \mu} \left[ {\rm e}^{{-}a \mu} Z_a(u_+) Z_a(u_-) \right]. \end{equation}

Observe that (2.30)–(2.38) only reference ${\rm I}_{a,b}$, and its derivatives with respect to $a$ and $b$, evaluated at $a = b = 1$. We thus set $a=b$ in (4.10) and, noting that $Z_a(u_\pm ) \to 0$ as $\mu \to \infty$ since $\mathrm {Im}(u_\pm ) > 0$, we find

(4.11)\begin{equation} {\rm I}_{a,a}^+= - \frac{1}{2 \sqrt{a} \zeta_d } Z_a(\zeta_+^+) Z_a(\zeta_-^+) , \end{equation}

where we have introduced

(4.12)\begin{equation} \zeta_\pm^+= \left. u_\pm \right|_{\mu = 0} = \frac{\sqrt[+]{1 + 8\zeta_d\zeta} \pm 1}{4\zeta_d}. \end{equation}

From (4.11) and (4.12), it is clear that ${\rm I}_{a,a}$ is a multivalued function with a branch point at $\zeta = -1/8\zeta _d$. Its second branch can be obtained by considering the $\sqrt [-]{}$ branch of the square root in (4.12). This means that both branches can be summarised by defining

(4.13)\begin{equation} \zeta_\pm^\lambda \equiv \frac{\sqrt[\lambda]{1 + 8\zeta_d\zeta} \pm 1}{4\zeta_d} = \frac{\lambda\sqrt[+]{1 + 8\zeta_d\zeta} \pm 1}{4\zeta_d}, \end{equation}

where $\lambda = \pm$ labels the branch. Therefore, ${\rm I}_{a,a}$ can be written as

(4.14)\begin{equation} {\rm I}_{a,a}^{\lambda} ={-} \frac{1}{2 \sqrt{a} \zeta_d } Z_a(\zeta_+^\lambda) Z_a(\zeta_-^\lambda). \end{equation}

Equation (4.14) is the key result of this paper. We shall henceforth refer to it as the generalised plasma dispersion function, in that it is the generalisation of the usual plasma dispersion function (3.1) to include the resonances associated with the magnetic drifts arising in a non-uniform magnetic field. In § 5.1, we show that, in the appropriate limits, the generalised plasma dispersion function reduces to the already-known solutions discussed in § 3. It is worth stressing that (4.14) is an exact result: no approximations have been made in deriving it from (2.55). Furthermore, the fact that (4.14) is composed of a product of plasma dispersion functions, for the evaluation of which there are numerous efficient algorithms, means that it is very fast to evaluate numerically. In § 5.2, we compare our expression for ${\rm I}_{a,b}$ with the numerical solver of Gürcan (Reference Gürcan2014).

It can be shown, via a similar procedure to that used to obtain (4.10) (see appendix A), that the related integral ${\rm J}_{a,b}^+$ (2.56) can be expressed exactly in terms of ${\rm I}_{a,b}^{+}$ and plasma dispersion functions as

(4.15)\begin{equation} \left( 1 - \frac{2b}{a} \right) {\rm J}_{a,b}^+= \frac{1}{2 a \zeta_d }\left[Z_a(\zeta_+^+) - Z_a(\zeta_-^+) \right] + \frac{b}{2a \zeta_d} {\rm I}_{a,b}^+, \end{equation}

and so

(4.16)\begin{equation} {\rm J}_{a, a}^\lambda ={-}\frac{1}{2a\zeta_d} \left[Z_a(\zeta_+^\lambda) - Z_a(\zeta_-^\lambda) \right] - \frac{1}{2\zeta_d} {\rm I}_{a, a}^\lambda. \end{equation}

It is crucial to realise that the $\lambda =+$ branch of the functions ${\rm I}_{a,a}^{\lambda }$ and ${\rm J}_{a,a}^{\lambda }$ is the ‘more important’ one, in the sense that it is the branch that is equal to the integrals (2.55) and (2.56) for $\mathrm {Im}(\zeta ) > 0$. Thus, it is also the branch that is used in the inverse Laplace transform over ${C_\sigma }$, as in (2.24). Therefore, we shall refer to the $\lambda =+1$ branch as the ‘principal’ branch of ${\rm I}_{a,a}^{\lambda }$ and ${\rm J}_{a,a}^{\lambda }$.

4.3 Derivatives of the generalised plasma dispersion function

In addition to (2.55) and (2.56), the matrix elements (2.30)–(2.38) require the partial derivatives of these expressions with respect to $a$ and $b$. There are two factors that conspire to simplify the necessary calculations. First, we only need ${\rm I}_{a,b}$, ${\rm J}_{a,b}$ and their derivatives at $a=b=1$. Second, the derivatives $\partial _a$ and $\partial _b$ often appear in the combination $\partial _a + \partial _b$. Notice that, by the chain rule,

(4.17)\begin{equation} (\partial_a + \partial_b)f_{a,b}\left\vert_{a=b=1} = \partial_a f_{a,a}\right\vert_{a=1} \end{equation}

for any (appropriately smooth) function $f$. Using this, we can rewrite (2.30)–(2.38) in a way that involves only ${\rm I}_{a, a}$, ${\rm J}_{a, a}$, $\partial _a {\rm I}_{a,b}\vert _{a=b}$, $\partial _b {\rm I}_{a,b}\vert _{a=b}$, $\partial _b^2 {\rm I}_{a,b}\vert _{a=b}$, and $\partial _b {\rm J}_{a,b}\vert _{a=b}$. For example,

(4.18)\begin{equation} {\mathsf{L}}_{\phi B} = \sum_s \frac{q_s n_{0s}}{q_r n_{0r}} \left[\zeta_s - \zeta_{*s} + \eta_s \zeta_{*s} \left( \partial_a + \frac{3}{2} \right) \right] \left.\left(\partial_b \left.{\rm I}_{a,b}^{(s)}\right|_{a=b}\right)\right|_{a=1}, \end{equation}

where we have also taken advantage of (2.51). Due to their unwieldy length, the calculations of the required derivatives of ${\rm I}_{a,b}$ and ${\rm J}_{a,b}$ are relegated to appendix B.

4.4 Branches of the dispersion function

Our choice of the principal branch cut for the square root gives the branches of the dispersion function $D$ (see § 2.2) several nice properties stemming from the relationship $\sqrt [+]{z^*} = \sqrt [+]{z}^*$ for any complex $z$. In appendix C, we show that ${\rm I}_{a, a}$ and ${\rm J}_{a, a}$ satisfy

(4.19)\begin{align} {\rm I}_{a, a}^\lambda(-\zeta^*, -\zeta_d) &={-}{\rm I}_{a, a}^\lambda(\zeta, \zeta_d)^*, \end{align}

(4.20)\begin{align}{\rm J}_{a, a}^\lambda(-\zeta^*, -\zeta_d) &= {\rm J}_{a, a}^\lambda(\zeta, \zeta_d)^*, \end{align}

and

(4.21)\begin{gather} {\rm I}_{a, a}^\lambda(\zeta^*, \zeta_d) = {\rm I}_{a, a}^{-\lambda}(\zeta, \zeta_d)^*, \end{gather}

(4.22)\begin{gather}{\rm J}_{a, a}^\lambda(\zeta^*, \zeta_d) = {\rm J}_{a, a}^{-\lambda}(\zeta, \zeta_d)^*. \end{gather}

Relations (4.19)–(4.22) are also valid for the $a$ and $b$ derivatives of ${\rm I}_{a,b}$ and ${\rm J}_{a,b}$. Of course, the functions ${\rm I}_{a,b}$ and ${\rm J}_{a,b}$ are double-valued for each of the species $s$, and so the dispersion function $D$ has $2^{N}$ branches for a system with $N$ species. Letting $\boldsymbol {\lambda } \equiv (\lambda _1, \lambda _2,\ldots, \lambda _N)$ be the vector of choices of the branch for each species, we can prove that $D$ satisfies (see appendix C)

(4.23)\begin{align} D^{\boldsymbol{\lambda}}(p^*, -\boldsymbol{k}) &= D^{\boldsymbol{\lambda}}(p, \boldsymbol{k})^* , \end{align}

(4.24)\begin{align}D^{\boldsymbol{\lambda}}({-}p^*, \boldsymbol{k}) &= D^{-{\boldsymbol{\lambda}}}(p, \boldsymbol{k})^*. \end{align}

These imply two different pairings of roots of the dispersion relation (2.48); see figures 11 and 12 in appendix C for a visual illustration of (4.23) and (4.24). Note that when using the superscript $\boldsymbol {\lambda }$, we are referring to a particular branch, while without it, $D$ refers to all branches simultaneously.

Equation (4.23) implies that solutions to the dispersion relation, i.e. $D=0$, come in pairs $(p, k_y) \leftrightarrow (p^*, -k_y)$, which is the condition for the fields $\phi$, $A_\parallel$, and $\delta B_\parallel$ to remain real for all $t$. Therefore, such a pairing is bound to exist for all roots of $D = 0$, i.e. when all branches are considered. The choice of the principal branch of the square root makes this pairing also valid within each individual branch of $D$, hence justifying our adoption of it in § 4.1. In § 6, we shall see that there is a better choice of branch for the purposes of performing the inverse Laplace transform.

Additionally, (4.24) says that if $p$ is a solution to $D=0$ for a given poloidal wavenumber $k_y$, then so is $-p^*$ for the same $k_y$ but for a different branch. At first glance, this might seem to imply that solutions to (2.48) always come in pairs, one stable and one unstable. While this is true if all branches of $D$ are considered, the time evolution given by the inverse Laplace transform (2.24) does not necessarily pick up contributions from all solutions to $D=0$; one cannot mix-and-match roots from different branches at will. In § 6, we shall see that the roots of (2.48) picked up by (2.24) depend on the choice of branch cut. However, only the principal branch, given by $\boldsymbol {\lambda } = (+, +,\ldots, +)$, contributes linearly unstable solutions.

5.1 Asymptotic expansions of ${\rm I}_{a,b}$

Let us now show that (4.14) asymptotes to the known limits discussed in § 3, as it should. First, in the limit of $\zeta _d \to 0$, i.e. the limit of vanishing magnetic curvature, we employ the expansions

(5.1)\begin{align} \zeta_+^+&= \frac{1}{2\zeta_d} + \zeta - 2\zeta_d\zeta^2 + {O}{\left(\zeta_d^2\right)} , \end{align}

(5.2)\begin{align}\zeta_-^+&= \zeta - 2\zeta_d\zeta^2 + 8\zeta_d^2\zeta^3 + {O}{\left(\zeta_d^3\right)} , \end{align}

and the asymptotic form $Z(\zeta ) \sim -\zeta ^{-1}$ for finite $\mathrm {Im}(\zeta )$ but $|\mathrm {Re}(\zeta )| \to \infty$, to find

(5.3)\begin{equation} Z_a(\zeta_+^+) \sim \frac{-2\zeta_d}{\sqrt{a}},\quad Z_a(\zeta_-^+) \sim Z_a(\zeta). \end{equation}

Therefore, in the limit $\zeta _d \to 0$, the principal branch satisfies

(5.4)\begin{equation} {\rm I}_{a,a}^+\sim\frac{1}{a}Z_a(\zeta), \end{equation}

in agreement with (3.2). This is visualised in figure 3. One can perform an analogous calculation with ${\rm J}_{a,b}^+$ to obtain the second expression in (3.2). Note that the second branch satisfies

(5.5)\begin{equation} {\rm I}_{a,a}^-\sim{-}\frac{1}{a}Z_a(-\zeta) \end{equation}

in the limit $\zeta _d \to 0$, which is not related to the correct expression for ${\rm I}_{a, a}$ at zero magnetic curvature. The ‘connection’ between the two branches, viz. the branch cut, is ‘sent to infinity’ as $\zeta _d \to 0$ (see figure 3), and so the second branch ${\rm I}_{a, a}^-$ is ‘lost’ in the limit of zero magnetic curvature. In this way, the dispersion function loses all but one of its branches and becomes single-valued.

Figure 3.

A plot of the principal branch ${\rm I}^+_{1,1}(\zeta, \zeta _d)$ in the complex plane for decreasing values of $\zeta _d$ (from left to right). The black cross denotes the branch point $\zeta = -1/8\zeta _d$. Panel (d) shows ${\rm I}_{1,1}^+(\zeta, 0) = Z(\zeta )$. As $\zeta _d \to 0^+$, the branch point, alongside the entire branch cut, is pushed towards $\mathrm {Re}(\zeta ) \to -\infty$. If $\zeta _d$ were negative, the branch cut would instead join the branch point with $\mathrm {Re}(\zeta )\to +\infty$, to which the branch cut would be pushed in the limit of $\zeta _d\to 0^-$.

Similarly, the 2-D limit can be found by taking the limit $\zeta \sim \zeta _d\to \infty$, which is equivalent to dropping the $u$ term from the denominators in (2.55) and (2.56). In this case,

(5.6)\begin{equation} \zeta_\pm^+= \sqrt[+]{\frac{\zeta}{2\zeta_d}} \pm \frac{1}{4\zeta_d} + O(\zeta_d^{{-}2}), \end{equation}

and so one obtains (3.4).

Figure 4 compares the exact expressions (4.14) and (4.16) with their known asymptotic limits in the case of vanishing magnetic drifts (3.2) and 2-D perturbations (3.4), respectively. It is evident that, while these known asymptotic limits are obtained in the cases of small and large $\zeta _d$, they are not a good approximation of ${\rm I}_{a,b}$ and ${\rm J}_{a,b}$ for $\zeta _d\sim 1$, as one would expect.

Figure 4.

Plots of the asymptotic convergence of ${\rm I}_{a,b}$ and ${\rm J}_{a,b}$ to their known limits. Panels (a, b) demonstrate the convergence of ${\rm I}_{1,1}^+$ and ${\rm J}_{1,1}^+$, given by (4.14) and (4.16), to their small- and large-$\zeta _d$ limits, given by (3.2) and (3.4), respectively. We define the relative difference of two functions $f$ and $g$ as $|f-g|/\min \{|f|, |g|\}$. Note that this is ill-defined if one of the functions is identically zero, so for the ${\rm J}_{1,1}$ comparison in panel (b), we plot simply $|{\rm J}_{1,1}|$ because we expect to recover ${\rm J}_{1,1} = 0$ in the 2-D limit (3.4). The solid and dotted lines in panels (a, b) show the average and maximum relative difference, respectively, as computed over a grid of $32\times 32$ points for $\zeta$, equally spaced in $\mathrm {Re}(\zeta ) \in [1, 1]$, $\mathrm {Im}(\zeta ) \in [0, 1]$, for each value of $\zeta _d$. Panel (c) demonstrates the convergence of the real (dashed) and imaginary (dash–dotted) parts of ${\rm I}_{1,1}^+$ to the small- and large-$\zeta _d$ limits for a fixed $\zeta = 1+{\rm i}$, which are given by $Z(\zeta )$ and $-Z^2(\sqrt {\varOmega })/\zeta _d$, with $\varOmega = \zeta /2\zeta _d$, respectively.

5.2 Numerical comparison with Gürcan (Reference Gürcan2014)

Gürcan (Reference Gürcan2014) consider a very similar problem to that on which this paper has focused but from a numerical perspective. In particular, they discuss the numerical integration of the function

(5.7)\begin{equation} \mathcal{I}_{nm}(\zeta_\alpha, \zeta_\beta, b) = \frac{2}{\sqrt{\rm \pi}}\int_0^\infty \,\mathrm{d} x_\perp\int_{-\infty}^\infty\,\mathrm{d} x_\parallel \: \frac{x_\perp^n x_\parallel^m {\rm J}_0^2(\sqrt{2b}x_\perp) {\rm e}^{{-}x_\parallel^2-x_\perp^2}}{x_\parallel^2 + x_\perp^2/2 + \zeta_\alpha - \zeta_\beta x_\parallel}, \end{equation}

defined for $\mathrm {Im}(\zeta _\alpha ) > 0$, real $\zeta _\beta$ and $b$, and $n>1$. With a few algebraic manipulations, it can be shown that, for odd $n$,

(5.8)\begin{equation} \mathcal{I}_{nm}\left(-\frac{\zeta}{2\zeta_d}, -\frac{1}{2\zeta_d}, 0\right) ={-}\frac{1}{\zeta_d}\left\{ \begin{array}{@{}ll} (-\partial_a)^{{m}/{2}}(-\partial_b)^{({n-1})/{2}} {\rm I}_{a,b}^+(\zeta, \zeta_d)\vert_{a=b=1} & \text{for }m\text{ even},\\ (-\partial_a)^{({m-1})/{2}}(-\partial_b)^{({n-1})/{2}} {\rm J}_{a,b}^+(\zeta, \zeta_d)\vert_{a=b=1} & \text{for }m\text{ odd}. \end{array} \right. \end{equation}

The above expression is actually correct only for $\zeta _d < 0$, otherwise (5.7) computes the second ($\lambda = -$) branch of ${\rm I}_{a,b}$, ${\rm J}_{a,b}$ and their derivatives. In the case $\zeta _d < 0$, the requirement $\mathrm {Im}(\zeta _\alpha ) > 0$ implies that $\mathrm {Im}(\zeta ) > 0$. Figure 5 shows a comparison between the values obtained via the results of this work [represented by (4.14), (4.16), (B12)–(B15)] and the Gürcan (Reference Gürcan2014) result (5.7) in the region $\mathrm {Re}(\zeta )\in [-10, 10]$, $\mathrm {Im}(\zeta )\in [0, 10]$, ${\zeta _d \in [-10, -0.001]}$. There is good agreement for all tested values of $\zeta$ and $\zeta _d$, with less than $1\,\%$ relative difference in most cases. It is important to stress that as our solution uses only standard functions, e.g. the plasma dispersion function $Z$ and $\sqrt {}$, for which there exist very efficient numerical algorithms. We found that even a naïve, unoptimised Python implementation took anywhere between $20$ and $80$ times less time to compute ${\rm I}_{a,b}$, ${\rm J}_{a,b}$ and their derivatives than the direct numerical integration of (5.7) implemented in Fortran at https://github.com/gurcani/zpdgen.

Figure 5.

Mean (solid) and maximum (dotted) relative difference (defined as in figure 4) between expressions (4.14), (4.16), (B12)–(B15) and their equivalents derived from (5.8), computed via the code published at https://github.com/gurcani/zpdgen. For each $\zeta _d$, we evaluated the respective functions at an equally spaced grid of $32\times 32$ points in the region ${\mathrm {Re}(\zeta )\in [-10, 10]}$, $\mathrm {Im}(\zeta )\in [0, 10]$.

7.1 Bessel functions

The drift-kinetic assumption is perhaps the more egregious approximation, especially given that the presence of finite-Larmor-radius effects can have a non-trivial impact on the plasma dynamics (see, e.g. Smolyakov et al. Reference Smolyakov, Yagi and Kishimoto2002; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020, Reference Parisi, Parra, Roach, Hardman, Schekochihin, Abel, Aiba, Ball, Barnes and Chapman-Oplopoiou2022, and references therein). Thankfully, however, it can be relaxed if one is willing to pay the price of complicated analytical expressions. Noting that $2 {\rm J}_0(b_s) {\rm J}_1(b_s) = - \partial {\rm J}_0^2(b_s)/\partial b_s$, it is clear that the Bessel functions ${\rm J}_0$ and ${\rm J}_1$ always appear quadratically in (2.39)–(2.43), for which there are known, rapidly converging Taylor series (Neumann Reference Neumann1871; Watson Reference Watson1966):

(7.1)\begin{equation} {\rm J}_n^2(b_s) = \sum_{m=0}^\infty \frac{({-}1)^m (2n+2m)!}{m!(2n+m)![(n+m)!]^2}\left(\frac{b_s}{2}\right)^{2n+2m}. \end{equation}

Using this expansion in (2.39)–(2.43), one can, in principle, compute each of the resulting integrals analytically and thus obtain an absolutely convergent series for the resulting gyrokinetic dispersion relation. This is done by noticing that their argument $b_s$ only appears quadratically as $b_s^2 = \mu k_\perp ^2 \rho _s^2$ and thus the additional factors of $\mu$ can be handled by partial differentiation with respect to $b$ before setting $a = b$ in (2.39)–(2.43). For example, (2.39) would give

(7.2)\begin{equation} \mathcal{I}_{a,b}^{(s)}(\zeta_s, \zeta_{ds}, \zeta_{ds}) = \sum_{m=0}^\infty \frac{(2m)!}{(m!)^4}\left(\frac{k_\perp \rho_s}{2}\right)^{2m} \partial_b^m {\rm I}_{a,b}(\zeta_s, \zeta_{ds}), \end{equation}

with the other required integrals, viz. $\mathcal {J}^{(s)}_{a,b}$ and $\partial _a \mathcal {I}^{(s)}_{a,b}$, satisfying similar expressions. We remind the reader that $\mathcal {I}_{a,b}^{(s)}$ refers to the FLR-containing integral (2.39), while ${\rm I}_{a,b}$ is the integral (2.55) on which we have focused throughout most of this paper. Doing this calculation by hand seems rather daunting given the complicated expressions even for the low-order derivatives $\partial _b^2 {\rm I}_{a,b}$ and $\partial _b {\rm J}_{a,b}$ [see (B14) and (B15), respectively]. In practice, however, only a few terms would be needed due to the rapid convergence of the Taylor series (7.1). Those wishing to compute these terms to an arbitrary order may want to do so by using symbolic libraries (e.g. those in Wolfram Mathematica) to calculate the derivatives analytically, which can then be imported into an associated numerical solver. An alternative approach would be to implement a recursive scheme to calculate numerically the $m$th-order derivatives from the $(m-1)$th ones.

7.2 General magnetic drifts

Our second approximation was to neglect the difference between the curvature and ${\boldsymbol {\nabla }} \! B$ drifts, taking their associated drift frequencies to be equal, i.e. $\zeta _{\kappa s} = \zeta _{{\boldsymbol {\nabla }} \! B s}$, as in (2.54). While this approximation is relatively well satisfied in the context of magnetic-confinement fusion, there are certainly other systems in which it is not, e.g. space and astrophysical plasmas. By a simple change of variables to $\mu ' = \zeta _{Bs}\mu /\zeta _{\kappa s}$ in (2.52), we find

(7.3)\begin{align} {\rm I}_{a,b}^{(s)}(\zeta_s, \zeta_{\kappa s}, \zeta_{Bs}) & = \frac{1}{\sqrt{\rm \pi}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu\: \frac{{\rm e}^{{-}a u^2 - b \mu}}{u - \zeta_s + \left(2u^2 \zeta_{\kappa s} + \mu \zeta_{Bs} \right)} \nonumber\\ & = \frac{\zeta_{\kappa s}}{\sqrt{\rm \pi}\zeta_{Bs}} \int_{-\infty}^\infty \,\mathrm{d} u \int_0^\infty \,\mathrm{d} \mu' \: \frac{{\rm e}^{{-a u^2 - (b\zeta_{\kappa s}/\zeta_{Bs})\mu'}}}{u - \zeta_s + \zeta_{\kappa s}\left(2u^2 + \mu' \right)}. \end{align}

Therefore, the integral (2.52) that enters (2.30)–(2.38) at $a=b$ can be found in terms of the known integrals (2.55) via

(7.4)\begin{equation} {\rm I}_{a,a}^{(s)}(\zeta_s, \zeta_{\kappa s}, \zeta_{Bs}) = \frac{\zeta_{\kappa s}}{\zeta_{Bs}} {\rm I}_{a,b}(\zeta_s, \zeta_{\kappa s}), \end{equation}

where now $b = a \zeta _{\kappa s}/\zeta _{Bs}$. Finally, to find ${\rm I}_{a,b}$ for $a\neq b$, one can Taylor expand

(7.5)\begin{equation} {\rm I}_{a,b} = \sum_{m=0}^\infty \frac{(b-a)^m}{m!} \left.\partial_b^m {\rm I}_{a,b}\right|_{a=b}. \end{equation}

Fortunately, just as in § 7.1, one can find closed, albeit complicated, analytical expressions for $\partial _b^m {\rm I}_{a,b}|_{a=b}$ for any $m$. The expansion should converge for arbitrary positive $a$ and $b$ since (7.5) is equivalent to expanding ${\rm e}^{-(b-a)\mu }$ in (4.10) using its absolutely convergent Taylor series.