An analytical form of the dispersion function for local linear gyrokinetics in a curved magnetic field

Abstract

Starting from the equations of collisionless linear gyrokinetics for magnetised plasmas with an imposed inhomogeneous magnetic field, we present the first known analytical, closed-form solution for the resulting velocity-space integrals in the presence of resonances due to both parallel streaming and constant magnetic drifts. These integrals are written in terms of the well-known plasma dispersion function (Faddeeva & Terent'ev, Tables of Values of the Function $w(z)=\exp (-z^2)(1+2i/\sqrt {\pi }\int _0^z \exp (t^2) \,\mathrm {d} t)$ for Complex Argument, 1954. Gostekhizdat. English translation: Pergamon Press, 1961; Fried & Conte, The Plasma Dispersion Function, 1961. Academic Press), rendering the subsequent expressions simpler to treat analytically and more efficient to compute numerically. We demonstrate that our results converge to the well-known ones in the straight-magnetic-field and two-dimensional limits, and show good agreement with the numerical solver by Gürcan (J. Comput. Phys., vol. 269, 2014, p. 156). By way of example, we calculate the exact dispersion relation for a simple electrostatic, ion-temperature-gradient-driven instability, and compare it with approximate kinetic and fluid models.

1 Introduction

The investigation of the linear-stability properties of magnetically confined plasmas is crucial for the design of magnetic-confinement-fusion devices. The heat and particle losses in these devices are dominated by turbulent fluctuations, which are themselves excited by linear instabilities driven by the gradients of the plasma equilibrium (Rudakov & Sagdeev 1961; Coppi, Rosenbluth & Sagdeev 1967; Pogutse 1968; Guzdar et al. 1983; Hugill 1983; Liewer 1985; Waltz 1988; Wootton et al. 1990; Cowley, Kulsrud & Sudan 1991; Kotschenreuther et al. 1995; Xanthopoulos et al. 2007; Ongena et al. 2016). In most cases, the strong toroidal magnetic field constrains the plasma fluctuations to have typical temporal scales that are slow compared to the frequency of the Larmor motion of the particles, and to be anisotropic in space: length scales along the magnetic field are comparable to the size of the device, while those perpendicular to it are comparable to the Larmor radii of the particles. Therefore, the plasma dynamics can often be treated using the gyrokinetic formalism (Frieman & Chen 1982; Sugama et al. 1996; Howes et al. 2006; Abel et al. 2013; Catto 2019).

When solving the linear gyrokinetic equation, one inevitably encounters resonant velocity-space integrals that need to be evaluated, analytically or numerically, to obtain the dispersion relation for the linear modes present within the system. The most basic of these resonances results from the parallel (to the magnetic field) streaming of particles, first discussed by Landau (1946). However, in the presence of an inhomogeneous equilibrium magnetic field, one is presented with a qualitatively different type of resonance due to the magnetic drifts of the particles. Evaluating these resonant integrals analytically, in the presence of both parallel streaming and magnetic drifts that are constant along the magnetic field, and without further approximations (such as those used in, e.g. Terry, Anderson & Horton 1982; Kim et al. 1994) has remained an open research question, despite some progress being made numerically (Gürcan 2014; Gültekin & Gürcan 2018; Gültekin & Gürcan 2020; Parisi et al. 2020). However, it is well known that there are instabilities that exist only in the presence of curved magnetic fields, e.g. the toroidal ion-temperature-gradient (ITG) instability (Pogutse 1968; Guzdar et al. 1983; Waltz 1988; Kotschenreuther et al. 1995). Often, such instabilities are the dominant ones in toroidal plasmas. Thus, the exact inclusion of the magnetic-drift resonance in the analytical theory of linear gyrokinetics is expected to lead to qualitative changes in the behaviour of the resulting dispersion relation and to allow for a more complete treatment of the linear-stability properties of strongly magnetised plasmas.

In this work, we present closed forms for the aforementioned resonant integrals. These are written in terms of the plasma dispersion function (Faddeeva & Terent'ev 1954; Fried & Conte 1961). They allow us to find a closed expression for the drift-kinetic dispersion relation, or an absolutely convergent series for the gyrokinetic one via Taylor expansions. The inclusion of magnetic drifts in the linear gyrokinetic problem introduces two distinct changes: (i) quantitatively, in that it significantly modifies the growth rates and frequencies of linear solutions; and (ii) qualitatively, by introducing a multivalued dispersion function. The latter has important consequences for the form of the dispersion relation and its solution, some of which have already been described in the literature (Kuroda et al. 1998; Sugama 1999).

The rest of the paper is organised as follows. We begin by summarising how the gyrokinetic dispersion relation, and the resonant velocity-space integrals of which it is comprised, emerge from the Fourier–Laplace transform of the linear gyrokinetic equations in § 2. Then, in § 3, we discuss already-known solutions for these integrals and the asymptotic limits in which they apply. The main result of this work is presented in § 4, where we derive the exact solution to one particular resonant integral – the ‘generalised plasma dispersion function’ – to which all others will be related. In § 5.1, we show both analytically and numerically that the generalised plasma dispersion function asymptotes to the known solutions in the cases of zero magnetic curvature and of two-dimensional (2-D) perturbations, while § 5.2 demonstrates that our expressions are in agreement with the numerical solver published by Gürcan (2014). Section 6 discusses the analytic continuation of these functions and the subsequent solution to the inverse-Laplace-transform problem by which we obtain the solution to the linear gyrokinetic system. In § 7, we show how the results obtained in § 4 can be generalised to the gyrokinetic case via absolutely convergent Taylor expansions. In § 8, we give an example calculation for the electrostatic ITG instability and compare it with known kinetic and fluid limits. Finally, our results are summarised and possible extensions discussed in § 9.

2 Collisionless gyrokinetic linear theory

In this section, we demonstrate how the resonant kinetic integrals that are the main focus of this paper emerge naturally from considerations of linear, collisionless local gyrokinetic theory with constant geometric coefficients [see the discussion following (2.10)]. Readers already familiar with gyrokinetic theory may wish to skip ahead to § 2.3, working backwards where further clarification is required.

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. 2013 or Catto 2019); 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 1995): 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. 2006; Ivanov et al. 2020; Adkins et al. 2022; Ivanov, Schekochihin & Dorland 2022) 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. 2013], 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. 2006)

(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 1972) 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 (1946), 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 rate¹

(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 Previous solutions

Before tackling the task of analytically integrating (2.55) and (2.56), we shall briefly discuss some special cases in which these expressions are already known within the literature. A reader already familiar with these solutions may wish to skip ahead to § 4, working backwards if further clarification is required.

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 1954; Fried & Conte 1961):

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

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 1971; Lee et al. 1987) and ion-temperature-gradient (see, e.g. Rudakov & Sagdeev 1961; Coppi et al. 1966; Sauter, Vaclavik & Skiff 1990; Brunner & Vaclavik 1998; Smolyakov, Yagi & Kishimoto 2002) 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. 2006).

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 1989), 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. 1984; Biglari et al. 1989; Kuroda et al. 1998; Sugama 1999; Ricci et al. 2006; Helander et al. 2011; Mishchenko, Plunk & Helander 2018; Zocco et al. 2018).

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 1996; Gürcan 2014; Gültekin & Gürcan 2018; Gültekin & Gürcan 2020; Parisi et al. 2020 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 1996; Parisi et al. 2020). 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 The generalised plasma dispersion function

In this section, we detail the method by which (2.55) and (2.56) can be expressed in terms of the plasma dispersion function (3.1), making the resultant expressions simpler to treat both analytically and numerically. When solving the integrals, we will assume that $p$ remains within the region of analyticity $\mathrm {Re}(p) \geqslant \sigma > 0$, with $\sigma$ defined after (2.23). The analytic continuation will be performed only after obtaining expressions for (2.55) and (2.56) in terms of known functions. In the main text, we present the integration of (2.55); all other required expressions follow directly from this single integral and have been relegated to appendices A and B due to their complexity. The remainder of this section proceeds as follows. Section 4.1 discusses the multivalued nature of the integrand of (2.55) before evaluating the integral over $u$ in terms of plasma dispersion function (3.1), allowing us, in § 4.2, to obtain a closed form expression for (2.55) upon evaluating the remaining integral over $\mu$. In § 4.3, we discuss how the $\partial _a$ and $\partial _b$ derivatives of (2.55) and (2.56) can be obtained, with detailed calculations relegated to appendix B. Then, in § 4.4, we discuss some important properties of (2.55) and (2.56).

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

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 Comparison with known results

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 (2014)

Gürcan (2014) 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 (2014) 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]$.

6 Analytic continuation for the inverse Laplace transform

Together, the expressions (4.14) and (4.16) for ${\rm I}_{a,a}$ and ${\rm J}_{a,a}$, respectively, along with the derivatives (B12)–(B15), allow us to calculate $\boldsymbol{\mathsf{L}}$ and hence the Laplace-transformed fields (2.47). Recall that to determine the evolution of the system as a function of time, we need to compute the inverse Laplace transform

(6.1)\begin{equation} \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{equation}

where the contour of integration ${C_\sigma }$ is once again as in figure 1, and we remind the reader that $\hat {\boldsymbol {\chi }}_{\boldsymbol {k}}(p)$ is given by

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

where the vector of initial conditions $\boldsymbol {G}$ is given by (2.44)–(2.46).

The results of § 4.2 show that the entries of $\boldsymbol{\mathsf{L}}$ have branch points at $\zeta _s = -1/8\zeta _{ds}$, or equivalently, at $p = p_s$, where

(6.3)\begin{equation} p_s \equiv \frac{{\rm i} k_\parallel^2 v_{{\rm th}s}^2}{8 \omega_{ds}}, \end{equation}

but are otherwise free of poles since, apart from the square roots and the associated branch cuts, they are composed of entire functions. Recall that we have defined the branches of the dispersion function $D$ using the principal branch of the square root in (4.3). Therefore, the relevant branch $D^{\boldsymbol {\lambda }}$ that enters the inverse Laplace transform is the principal branch given by $\boldsymbol {\lambda } = (+,\ldots, +)$. This has branch cuts that connect the branch points $\zeta _s = -1/8\zeta _{ds}$ to $\zeta _s \to -\mathrm {sgn}(\zeta _{ds}) \infty$, or, equivalently, $p = p_s$ to $p \to {\rm i}\mathrm {sgn}(\zeta _{ds})\infty$. While this choice of the principal branch and branch cuts was convenient for obtaining the closed forms of ${\rm I}_{a,b}$, ${\rm J}_{a,b}$, and $D$, and their properties, it is not necessarily the best one for performing the inverse Laplace transform (6.1). Instead, we would like to rotate the branch cuts by $\mathrm {sgn}(\zeta _{ds}) {\rm \pi}/2$ around $p_s$, so that they are parallel to the real $p$ axis, as shown in figure 6. Let us call the branch $\mathcal {D}$ of the dispersion function obtained this way the ‘dispersion’ branch. Crucially, the rotation of the branch cuts does not disturb the values of the dispersion function at $\mathrm {Re}(p) > 0$. Therefore, $\mathcal {D}(p) = D^{(+,\ldots +)}(p)$ for $\mathrm {Re}(p) > 0$. This ensures that the ‘unphysical’ unstable zeros of the other branches of the dispersion function, which are a consequence of (4.24), do not contribute to the solution (see also discussion in § 4.4); the only unstable solutions that are picked up by the inverse Laplace transform are those of the principal branch.

Figure 6. This diagram shows the ‘principal’ (in blue) and ‘dispersion’ (in black) branch cuts for a plasma with one negatively and one positively charged species, labelled as $s_1$ and $s_2$, respectively.

With this choice for the branch cuts of the dispersion function, we are ready to perform the inverse Laplace transform (6.1). This is done in the usual way, viz. by pushing the integration contour ${C_\sigma }$ towards $\mathrm {Re}(p) \to -\infty$, with the proviso that it must be deformed so as not to cross any singularities, e.g. poles or branch cuts. Pushing the contour to the vertical line at $\mathrm {Re}(p) = \rho$, we find the new integration contour ${C_\rho }$ (see figure 7). Since there are no singularities between ${C_\sigma }$ and ${C_\rho }$, Cauchy's integral theorem ensures that the integrals over these two contours are equal. Taking the limit of $\rho \to -\infty$, it is evident that the contributions arising from the vertical segments of ${C_\rho }$ are exponentially small [they are exponentially small at any $t > 0$ because the integrand of the inverse Laplace transformation (2.24) contains a factor ${\rm e}^{\rho t}$], while those arising from the integration along the horizontal segments leading towards and away from the poles cancel, leaving the contributions from the poles. The integration around the branch cuts is more subtle and will be discussed shortly.

Figure 7. Same as in figure 1, except that the contour associated with the inverse Laplace transformation (6.1) has now been shifted to $\mathrm {Re}(p) = \rho$, deforming it such that it does not cross any of the poles or the branch cut. We denote this new contour ${C_\rho }$. The original contour is shown by the vertical dashed line. The integrals along ${C_\sigma }$ and ${C_\rho }$ are equal by Cauchy's integral theorem.

There are several singularities present in (6.2) and hence in the integrand in (6.1). The first is the so-called ‘ballistic response’ associated with the initial conditions contained within $\boldsymbol {G}$, arising from simple poles located along $\mathrm {Re}(p) = 0$, viz.

(6.4)\begin{equation} \lim_{\rho \rightarrow -\infty} \frac{1}{2 {\rm \pi}{\rm i}} \int_{{C_\rho}} \,\mathrm{d} p \: {\rm e}^{pt} \frac{{g}_{s \boldsymbol{k}}}{p + {\rm i} k_\parallel v_\parallel{+} {\rm i} \omega_{Ds}} = {g}_{s \boldsymbol{k}} {\rm e}^{{-{\rm i}(k_\parallel v_\parallel{+} \omega_{Ds})t}}, \end{equation}

where we have assumed that ${g}_{s \boldsymbol {k}}$ is a smooth function. Plugging this into (2.44)–(2.46), we find that the contribution to $\boldsymbol {\chi }_{\boldsymbol {k}}(t)$ due to the ballistic response can be written as

(6.5)\begin{gather} \boldsymbol{\chi}_{\boldsymbol{k}0}(t) = \sum_s \int \,\mathrm{d}^3 \boldsymbol{v} \ \boldsymbol{\mathsf{L}}^{{-}1}(-{\rm i}k_{\parallel} v_{\parallel}{-}{\rm i}\omega_{Ds}) {g}_{s \boldsymbol{k}} {\rm e}^{{-{\rm i}(k_{\parallel} v_{\parallel}{+} \omega_{Ds})t}} \nonumber\\ \left( \frac{q_s n_{0s}}{q_r n_{0r}} \frac{1}{n_{0s}} {\rm J}_0 (b_s), \frac{q_s n_{0s} v_{{\rm th}s}}{q_r n_{0r} v_{{\rm th}r} } \frac{1}{n_{0s}}\frac{v_\parallel}{v_{{\rm th}s}} {\rm J}_0 (b_s), - \frac{\beta_s}{2}\frac{1}{n_{0s}}\frac{v_\perp^2}{v_{{\rm th}s}^2} \frac{2 {\rm J}_1(b_s)}{b_s} \right)^{\rm T}. \end{gather}

There is a wealth of interesting physics that can arise from the ballistic response, see, e.g. Ewart et al. (2022) and references therein, in the context of the Vlasov–Poisson system. However, this is not the focus of the present work and so will not be discussed further.

Another source of non-analyticity are the solutions to the dispersion relation $\mathcal {D} = 0$, should any of these exist. The contributions to (6.1) arising from the zeros $p = p_j$ of $\mathcal {D}$ can be written as

(6.6)\begin{equation} \sum_j \mathrm{Res}[\hat{\boldsymbol{\chi}}_{\boldsymbol{k}}(p), p_j] {\rm e}^{p_j t}. \end{equation}

It is evident that unlike the ballistic response, whose time dependence is an oscillating exponential, the terms (6.6) can, in general, be exponentially decaying (i.e. stable) for $\mathrm {Re}(p_j) < 0$ or growing (i.e. unstable) for $\mathrm {Re}(p_j) > 0$.

Finally, singularities may arise from the functions (2.39)–(2.43) that are contained within both $\text {adj}\:\boldsymbol{\mathsf{L}}$ and $\det \boldsymbol{\mathsf{L}}$. As discussed above, these functions are free of poles, but are multivalued. Deforming the integration contour ${C_\rho }$ around their branch cuts (see figure 7) gives a non-trivial contribution to (6.1). Letting $\boldsymbol {B}_s(t)$ be the contribution from the integral around the branch cut due to a given species $s$, we can finally write the full solution for $\boldsymbol {\chi }_{\boldsymbol {k}}(t)$ as

(6.7)\begin{equation} \boldsymbol{\chi}_{\boldsymbol{k}}(t) = \boldsymbol{\chi}_{\boldsymbol{k}0}(t) + \sum_j \mathrm{Res}[\hat{\boldsymbol{\chi}}_{\boldsymbol{k}}(p), p_j] {\rm e}^{p_j t} + \sum_s \boldsymbol{B}_s(t). \end{equation}

In appendix D, we show that, in the long-time limit $t \to \infty$, the branch-cut contribution $\boldsymbol {B}_s$ for each species is dominated by that arising from the branch point itself, and exhibits an algebraic decay $\propto t^{-3/2}$. The same algebraic decay was found by Kim et al. (1994) and Kuroda et al. (1998) in their treatment of the toroidal ITG mode. Such a ‘continuum mode’ (Kuroda et al. 1998; Sugama 1999) is a direct consequence of the multivaluedness of (2.55) and (2.56), in that such multivaluedness gives rise to a branch point and to the resulting discontinuity. This behaviour is qualitatively different from that of a plasma in a straight magnetic field, whose dispersion function is single-valued, meaning that there are no branch cuts and hence no continuum modes. Note that non-exponentially decaying solutions to similar initial-value problems can also be found in other contexts; see, e.g. Taylor (1965) and Sedlàček (1995).

Equation (6.7) is our final expression for the time evolution of $\boldsymbol {\chi }_{\boldsymbol {k}}(t)$. Depending on whether there are any unstable solutions, we find that either: (i) there are solutions to $\mathcal {D}(p) = 0$ for $\mathrm {Re}(p) > 0$, in which case, the long-time solution is dominated by the solution with largest $\mathrm {Re}(p)$; or (ii) there are no solutions to $\mathcal {D}(p) = 0$ for $\mathrm {Re}(p) > 0$, in which case, the long-time solution is dominated by the ballistic response (6.4) and by waves with frequencies $\omega = {\rm i}p_s = -k_\parallel ^2v_{{\rm th}s}^2/8\omega _{ds}$ that exhibit a non-exponential decay $\propto t^{-3/2}$.

7 From drift kinetics to gyrokinetics

The analytical forms of the integrals derived in § 4 are not without their limitations: in their derivation, we assumed both the drift-kinetic limit and the case of equal magnetic drifts (see § 2.3). We will now devote some space to a brief discussion of how one can relax these assumptions.

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. 2002; Parisi et al. 2020, 2022, 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 1871; Watson 1966):

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

8 Electrostatic ITG: a detailed example

To illustrate the results of § 4, we provide an explicit calculation of the dispersion relation in the simple case of an electrostatic, ion-scale, temperature-gradient-driven instability, and compare the solution with well-known kinetic and fluid limits.

In particular, we consider a two-species plasma of ions and electrons of comparable temperatures, $T_{0i}\sim T_{0e}$. Since we want to consider electrostatic physics, we assume

(8.1)\begin{equation} \beta_s \sim (k_\perp d_s)^{{-}2} \ll 1, \end{equation}

where $d_s = \sqrt {m_sc^2 / 4{\rm \pi} n_{0s} q_s^2}$ is the skin depth of species s. Therefore, to lowest order, $A_{\parallel \boldsymbol {k}}$ and ${\delta \! B}_{\parallel \boldsymbol {k}}$ do not contribute to (2.17), and (2.29) simplifies to

(8.2)\begin{equation} {\mathsf{L}}_{\phi \phi} \frac{q_r {\hat{\phi}}_{\boldsymbol{k}}}{T_{0r}} + \text{G}_\phi = 0. \end{equation}

This implies that the dispersion relation is given simply by $\text {L}_{\phi \phi } = 0$. Furthermore, we consider the frequencies of the perturbations to be comparable to the parallel streaming and drift frequencies of the ions, as well as the magnetic-drift frequency, viz.

(8.3)\begin{equation} p\sim k_{\parallel} v_{{\rm th}i}\sim\omega_{di}\sim\omega_{*i}\sim\eta_i \omega_{*i}. \end{equation}

The relevant equilibrium length scales in our problem are thus the ion-density and ion-temperature gradients $L_{n_i}^{-1}$ and $L_{T_i}^{-1}$, respectively [see (2.6)], and the gradient of the magnetic field $L_B^{-1} \equiv -\partial \ln B_0 / \partial x$. In the small-mass-ratio limit, $m_e/m_i \ll 1$, (8.3) implies

(8.4)\begin{equation} k_\parallel v_{{\rm th}e} \sim \sqrt{\frac{m_i}{m_e}} k_\parallel v_{{\rm th}i} \gg p, \end{equation}

i.e. the electrons stream quickly along the fields lines. Thus, $\zeta _e \sim \zeta _{*e} \ll \zeta _i \sim \zeta _{*i}$, and the electron contributions to $\text {L}_{\phi \phi }$ can be ignored. Choosing $q_r = q_i = Ze$, $T_{0r} = T_{0i}$, and ${n_{0r} = n_{0i}}$, the expression (2.30) simplifies to

(8.5)\begin{equation} -{\mathsf{L}}_{\phi\phi} = 1 + \tau + \left[ \zeta_i - \zeta_{*i} + \eta_i \zeta_{*i} \left(\partial_a + \frac{3}{2}\right) \right] \left. {\rm I}_{a, a}^{(i)} \right|_{a=1}, \end{equation}

where $\tau \equiv T_{0i}/ZT_{0e}$ is the temperature ratio. To avoid carrying around an extra minus sign, we shall define $D \equiv -\text {L}_{\phi \phi }$, the object in whose zeros we shall be interested. Using (4.14), we obtain the principal branch of the ITG dispersion relation:

(8.6)\begin{equation} D = 1 + \tau - \frac{\zeta - \zeta_*}{2 \zeta_d} Z_+Z_-{+} \frac{\eta \zeta_*}{2\zeta_d} \left[ \left(\zeta_{+} Z_-{+} \zeta_- Z_+\right) + \left(\frac{\zeta}{\zeta_d} + \frac{1}{4\zeta_d^2} - 1\right) Z_+Z_- \right] = 0, \end{equation}

where we have dropped the $i$ subscripts, $\zeta _\pm$ are given by (4.12) and we are using the shorthand notation $Z_\pm \equiv Z(\zeta _\pm )$. Note that the principal branch (i.e. $\lambda =+$) is implicitly used everywhere, but we have dropped the associated superscripts to reduce the notational clutter.

We can use (5.2) and (5.6) to verify that (8.6) converges to the correct limits in the case of: vanishingly small magnetic gradients (i.e. $\zeta _d \to 0$)

(8.7)\begin{equation} D_\text{slab} = 1 + \tau + (\zeta - \zeta_*)Z(\zeta) - \eta \zeta_* \left[\zeta + \zeta^2Z(\zeta) - \frac{1}{2}Z(\zeta)\right] = 0; \tag{8.7}\end{equation}

and of 2-D perturbations (i.e. $\zeta \sim \zeta _*\sim \zeta _d\to \infty$)

(8.8)\begin{equation} D_\text{2D} = 1 + \tau - (\varOmega - \varOmega_*)Z(\sqrt{\varOmega})^2 + \eta\varOmega_*\left[2 \sqrt{\varOmega} Z(\sqrt{\varOmega}) + \left(2\varOmega - 1\right)Z(\sqrt{\varOmega})^2\right]= 0, \end{equation}

where $\varOmega = \zeta /2\zeta _d = ip/2\omega _d$ and $\varOmega _* = \zeta _* / 2\zeta _d = \omega _*/2\omega _d$. Note that (8.8) agrees with the expressions obtained by Biglari et al. (1989) and Zocco et al. (2018) in a similar limit to (8.3).

In figure 8, we compare the solutions to (8.6) and (8.7) for the case of zero density gradient, viz. $\omega _{*} = 0$, but non-zero temperature gradient, so $\eta \omega _{*} \propto L_{T_{i}}^{-1} \neq ~0$. The growth rates agree well only at simultaneously large perpendicular and small parallel wavelengths; this is to be expected given that the slab dispersion relation (8.7) does not capture the effect of magnetic drifts, which are most important at large parallel wavelengths. There is poorer agreement between the frequencies of the two dispersion relations.

Figure 8. A comparison between the growth rate and frequency of the most unstable solution to the kinetic dispersion relation with magnetic effects (8.6) and the slab dispersion relation (8.7), represented by the solid and dotted lines, respectively. Here, $\rho _s = \rho _i / \sqrt {2\tau }$ is the ion sound radius, and we have set $\tau =0.1$ and $\tau L_B / 2L_{T_{i}} = 2$.

We can also compare the solutions to (8.6) with those obtained from a simple three-field fluid model of the ITG instability in a slab with magnetic curvature. The model consists of the following equations:

(8.9)\begin{align}&\tau\frac{\partial {\varphi}}{\partial t} + \frac{\partial {u_\parallel}}{\partial z} - \frac{\rho_i v_{{\rm th}i}}{L_B} \frac{\partial {}}{\partial y}\left[(1+\tau)\varphi + \frac{\delta T_i}{T_{0i}}\right] = 0, \end{align}

(8.10)\begin{align}&\frac{\partial {u_\parallel}}{\partial t} + \frac{v_{{\rm th}i}^2}{2}\frac{\partial {}}{\partial z}\left[(1+\tau)\varphi + \frac{\delta T_i}{T_{0i}}\right] - \frac{2\rho_iv_{{\rm th}i}}{L_B} \frac{\partial {u_\parallel}}{\partial y} = 0, \end{align}

(8.11)\begin{align}&\frac{\partial {}}{\partial t}\frac{\delta T_i}{T_{0i}} + \frac{2}{3}\frac{\partial {u_\parallel}}{\partial z} + \frac{\rho_i v_{{\rm th}i}}{2 L_{T_i}} \frac{\partial {\varphi}}{\partial y} - \frac{2}{3}\frac{\rho_iv_{{\rm th}i}}{L_B} \frac{\partial {}}{\partial y}\left[(1+\tau)\varphi + \frac{7}{2}\frac{\delta T_i}{T_{0i}}\right] = 0, \end{align}

where $\varphi \equiv Ze\phi /T_{0i}$, $u_\parallel$, and $\delta T_i/T_{0i}$ are the perturbed electrostatic potential, ion parallel flow, and ion temperature, respectively. These equations can be derived by substituting a perturbed Maxwellian for $h_i$ in the ion gyrokinetic equation and taking the three relevant velocity moments (cf. Newton, Cowley & Loureiro 2010 or the cold-ion fluid model of Ivanov et al. 2022 but with additional $\tau \sim 1$ terms). Figure 9 shows a comparison between the kinetic and fluid growth rates at a fixed value of $\tau$ and varying $k_\parallel L_B$. We see that the fluid approximation is decent for small $k_\parallel L_B$, but fails for larger ones because of its lack of kinetic effects. Making the ions cold, i.e. lowering $\tau$, improves the accuracy of the fluid approximation, as in figure 10.

Figure 9. A comparison between the growth rate and frequency of the most unstable solution to the kinetic dispersion relation (8.6) and that obtained from the fluid equations (8.9)–(8.11), represented by the solid and dotted lines, respectively. The parameters used are the same as in figure 8.

Figure 10. A comparison between the growth rate and frequency of the most unstable solution to the kinetic dispersion relation (8.6) and that of the fluid equations (8.9)–(8.11), represented by the solid and dotted lines. Here, we have set $k_\parallel L_B = 1$ and $\tau L_B / 2L_{T_{i}} = 2$.

9 Summary and discussion

We have considered the problem of local linear gyrokinetics in a curved magnetic field, expressing the associated dispersion relation in terms of velocity-space integrals featuring resonances arising both from parallel streaming and from magnetic drifts (§ 2). Previously, exact solutions for these integrals were known either in the absence of magnetic drifts – leading to the well-known plasma dispersion function $Z (\zeta )$ – or in the 2-D limit (§ 3). In the case of drift kinetics (i.e. no finite-Larmor-radius effects) and equal magnetic drifts, we showed that these resonances can in fact be handled simultaneously without any additional approximations or expansions, and that the integrals can be expressed exactly in terms of a generalised plasma dispersion function consisting of products of $Z$ functions and their derivatives (§ 4). Since there exist known algorithms for the computation of the $Z$ function, the resulting expressions are efficient to evaluate numerically and can easily be handled analytically through known asymptotic expansions. Solutions to the exact dispersion relation for the electrostatic ITG instability, derived using this method, were then compared with approximate solutions in the previously known limits, showing poor agreement for the majority of parameters and wavenumbers considered (§ 8). This demonstrates that to properly capture the growth rate and frequency of kinetic instabilities in the presence of a curved magnetic field, one must simultaneously resolve the resonances associated with parallel streaming and magnetic drifts, for which this paper provides the first known exact analytical solution.

In § 7, we discussed how the assumptions of no finite-Larmor-radius effects and equal magnetic drifts can be relaxed using absolutely convergent Taylor-series expansions, and thus solve the more general linear gyrokinetic system. This results in expressions that naturally capture the multivaluedness of the underlying dispersion relation and handle the integration of resonant denominators exactly.

An immediate practical application of this work would be to use the derived analytical expressions to implement an efficient and accurate solver for drift-kinetic/gyrokinetic instabilities in the local limit considered in this paper. Such a solver could be used to benchmark both reduced models and gyrokinetic solvers. It could also be exploited to explore the equilibrium parameter space in search of new instabilities or to investigate the properties of subdominant ones, i.e. those whose growth rate is smaller than the largest growth rate in the system; this is typically difficult to do in most gyrokinetic solvers. Such subdominant instabilities have been proposed as one of the possible explanations for the lack of saturation observed in certain electromagnetic gyrokinetic simulations. With this in mind, we consider the implementation of such a gyrokinetic dispersion-relation solver to be a natural extension of this work that will produce a useful practical tool in the study of gyrokinetic instabilities and turbulence.