3.1 Integrating near poles

A challenge concerning the numerical integration is the treatment of the poles that occur in the term proportional to $\unicode[STIX]{x1D64F}_{n}$ in (2.9). The integrals in question are of the form

(3.1) $$\begin{eqnarray}\displaystyle I(p_{\bot })\equiv \int _{-\infty }^{+\infty }\text{d}p_{\Vert }\frac{\unicode[STIX]{x1D6FA}_{j}U\unicode[STIX]{x1D64F}_{n}}{\unicode[STIX]{x1D714}-k_{\Vert }v_{\Vert }-n\unicode[STIX]{x1D6FA}_{j}}\equiv \int _{-\infty }^{+\infty }\text{d}p_{\Vert }G(p_{\bot },p_{\Vert }) & & \displaystyle\end{eqnarray}$$

for $\unicode[STIX]{x1D6FE}>0$ . For sufficiently small $\unicode[STIX]{x1D6FE}$ , the denominator in (3.1) can become very small along the real $p_{\Vert }$ axis so that the grid sampling leads to large numerical errors in the integration. To describe how we evaluate these integrals, we first rewrite the integral in (3.1) in the more generic form

(3.2) $$\begin{eqnarray}\displaystyle {\mathcal{I}}\equiv \int _{-\infty }^{+\infty }\text{d}x\frac{g(x)}{x-t_{\text{r}}-\text{i}t_{\text{i}}}, & & \displaystyle\end{eqnarray}$$

where $x$ , $t_{\text{r}}$ and $t_{\text{i}}$ are real, $g(x)$ is a smooth function and the integration is performed along the real axis. We choose a symmetric interval $[t_{\text{r}}-\unicode[STIX]{x1D6E5},t_{\text{r}}+\unicode[STIX]{x1D6E5}]$ around $t_{\text{r}}$ where $\unicode[STIX]{x1D6E5}\ll g(t_{\text{r}})/g^{\prime }(t_{\text{r}})$ , and write

(3.3) $$\begin{eqnarray}\displaystyle {\mathcal{I}}=\int _{t_{\text{r}}-\unicode[STIX]{x1D6E5}}^{t_{\text{r}}-\unicode[STIX]{x1D6E5}}\text{d}x\frac{g(x)}{x-t_{\text{r}}-\text{i}t_{\text{i}}}+\text{rest}, & & \displaystyle\end{eqnarray}$$

where ‘rest’ refers to the integration outside the interval $[t_{\text{r}}-\unicode[STIX]{x1D6E5},t_{\text{r}}+\unicode[STIX]{x1D6E5}]$ . We define a function $f(x)$ to be odd with respect to $t_{\text{r}}$ if $f(x)=-f(2t_{\text{r}}-x)$ , and even with respect to $t_{\text{r}}$ if $f(x)=f(2t_{\text{r}}-x)$ . Following Longman (Reference Longman1958) and Davis & Rabinowitz (Reference Davis and Rabinowitz1984), we then separate the integrand into its odd and even parts with respect to $t_{\text{r}}$ as

(3.4) $$\begin{eqnarray}\displaystyle {\mathcal{I}} & = & \displaystyle \frac{1}{2}\int _{t_{\text{r}}-\unicode[STIX]{x1D6E5}}^{t_{\text{r}}+\unicode[STIX]{x1D6E5}}\text{d}x\left[\frac{g(x)}{x-t_{\text{r}}-\text{i}t_{\text{i}}}-\frac{g(2t_{\text{r}}-x)}{-x+t_{\text{r}}-\text{i}t_{\text{i}}}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2}\int _{t_{\text{r}}-\unicode[STIX]{x1D6E5}}^{t_{\text{r}}+\unicode[STIX]{x1D6E5}}\text{d}x\left[\frac{g(x)}{x-t_{\text{r}}-\text{i}t_{\text{i}}}+\frac{g(2t_{\text{r}}-x)}{-x+t_{\text{r}}-\text{i}t_{\text{i}}}\right]+\text{rest}.\end{eqnarray}$$

The integrand in the first integral in (3.4) is odd with respect to $t_{\text{r}}$ and thus vanishes after the integration over the symmetric interval around $t_{\text{r}}$ . The second integral, on the other hand, is even with respect to $t_{\text{r}}$ and thus

(3.5) $$\begin{eqnarray}\displaystyle {\mathcal{I}}=\int _{t_{\text{r}}}^{t_{\text{r}}+\unicode[STIX]{x1D6E5}}\text{d}x\left[\frac{g(x)}{x-t_{\text{r}}-\text{i}t_{\text{i}}}-\frac{g(2t_{\text{r}}-x)}{x-t_{\text{r}}+\text{i}t_{\text{i}}}\right]+\text{rest}. & & \displaystyle\end{eqnarray}$$

We define $\unicode[STIX]{x1D6E5}$ through a user-defined parameter $M_{\text{I}}$ so that $\unicode[STIX]{x1D6E5}\equiv M_{\text{I}}\,\unicode[STIX]{x0394}p_{\Vert }$ , where $\unicode[STIX]{x0394}p_{\Vert }$ is the size of a grid step in the parallel direction. We then define $\unicode[STIX]{x1D6FF}\equiv \unicode[STIX]{x1D6E5}/M_{\text{P}}$ , where $M_{\text{P}}$ is another user-defined parameter. Except for cases in which $|t_{\text{i}}|$ is extremely small, we apply a trapezoidal integration over $M_{\text{P}}$ steps of width $\unicode[STIX]{x1D6FF}$ to the integral in (3.5). The smoothness of $g(x)$ allows us to expand $g(x)$ around the nearest grid point of the $M_{\text{I}}$ grid points in the interval $[t_{\text{r}},t_{\text{r}}+\unicode[STIX]{x1D6E5}]$ using a Taylor series. By taking $\unicode[STIX]{x0394}p_{\Vert }$ to be sufficiently small, we can retain just the first two terms in the series without losing significant accuracy. Since the integral in (3.5) does not converge numerically if $|t_{\text{i}}|$ is extremely small, we implement the following procedure when $|t_{\text{i}}|\leqslant t_{\text{lim}}$ , where $t_{\text{lim}}$ is a user-defined parameter. We first rewrite (3.5) using truncated Taylor expansions of $g(x)$ and $g(2t_{\text{r}}-x)$ around $x=t_{\text{r}}$ as

(3.6) $$\begin{eqnarray}\displaystyle {\mathcal{I}}=\int _{t_{\text{r}}}^{t_{\text{r}}+\unicode[STIX]{x1D6E5}}\text{d}x\left[\frac{2\text{i}t_{\text{i}}g(t_{\text{r}})}{(x-t_{\text{r}})^{2}+t_{\text{i}}^{2}}+\frac{2g^{\prime }(t_{\text{r}})(x-t_{\text{r}})^{2}}{(x-t_{\text{r}})^{2}+t_{\text{i}}^{2}}\right]+\text{rest}. & & \displaystyle\end{eqnarray}$$

We determine $g(t_{\text{r}})$ and $g^{\prime }(t_{\text{r}})$ through linear interpolation between the neighbouring grid points to $t_{\text{r}}$ . The term proportional to $g^{\prime }(t_{\text{r}})$ in (3.6) converges numerically for any value of $t_{\text{i}}$ . We set the term proportional to $g(t_{\text{r}})$ equal to its small- $t_{\text{i}}$ limit, namely

(3.7) $$\begin{eqnarray}\displaystyle \int _{t_{\text{r}}}^{t_{\text{r}}+\unicode[STIX]{x1D6E5}}\text{d}x\frac{2\text{i}t_{\text{i}}g(t_{\text{r}})}{(x-t_{\text{r}})^{2}+t_{\text{i}}^{2}}=\text{i}\unicode[STIX]{x03C0}g(t_{\text{r}})\text{sgn}(t_{\text{i}}). & & \displaystyle\end{eqnarray}$$

We use this method for both the integration of $\unicode[STIX]{x1D74C}_{j}$ near poles and the principal-value integration that is necessary if $\unicode[STIX]{x1D6FE}=0$ .

3.2 Analytic continuation

If $\unicode[STIX]{x1D6FE}\leqslant 0$ , the integration in (2.9) requires an analytic continuation into the complex plane. If $f_{0j}$ were given as a closed algebraic expression, the analytic continuation would simply entail the evaluation of $f_{0j}(p_{\bot },p_{\Vert })$ at a complex value for $p_{\Vert }$ in the non-relativistic case. In our case, however, $f_{0j}$ is only defined on a real grid in $p_{\bot }$ and $p_{\Vert }$ , yet the analytic continuation of $f_{0j}$ is still uniquely defined. This leads to the known mathematical problem of numerical analytic continuation (Cannon & Miller Reference Cannon and Miller1965; Reichel Reference Reichel1986; Fujiwara et al. Reference Fujiwara, Imai, Takeuchi and Iso2007; Fu et al. Reference Fu, Zhang, Cheng and Ma2012; Zhang & Ma Reference Zhang and Ma2013; Kranich Reference Kranich2014). Our solution for this problem is our hybrid analytic continuation scheme. We note that this approach is only relevant for damped modes, i.e. $\unicode[STIX]{x1D6FE}\leqslant 0$ .

Landau’s rule of integration around singularities (Landau Reference Landau1946; Lifshitz & Pitaevskii Reference Lifshitz and Pitaevskii1981) leads to the following three cases with the appropriate residues for the evaluation of $I(p_{\bot })$ for general $\unicode[STIX]{x1D6FE}$ :

(3.8) $$\begin{eqnarray}\displaystyle I(p_{\bot })=\int _{C_{\text{L}}}\text{d}p_{\Vert }G(p_{\bot },p_{\Vert })=\left\{\begin{array}{@{}ll@{}}\displaystyle \int _{-\infty }^{+\infty }\text{d}p_{\Vert }G(p_{\bot },p_{\Vert })\quad & \text{if }\unicode[STIX]{x1D6FE}>0,\\ \displaystyle {\mathcal{P}}\int _{-\infty }^{+\infty }\text{d}p_{\Vert }G(p_{\bot },p_{\Vert })+\text{i}\unicode[STIX]{x03C0}\sum \text{Res}_{A}(G)\quad & \text{if }\unicode[STIX]{x1D6FE}=0,\\ \displaystyle \int _{-\infty }^{+\infty }\text{d}p_{\Vert }G(p_{\bot },p_{\Vert })+2\text{i}\unicode[STIX]{x03C0}\sum \text{Res}_{A}(G)\quad & \text{if }\unicode[STIX]{x1D6FE}<0,\end{array}\right. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $C_{\text{L}}$ is the contour of the Landau integration, which lies below the complex poles in the integrand. The integrations on the right-hand side of (3.8) are performed along the real axis, and ${\mathcal{P}}$ indicates the principal-value integral. The sum sign indicates the summation over the residues of all poles $A$ of the function $G$ . In a non-relativistic plasma, $G$ has one simple pole, and thus

(3.9) $$\begin{eqnarray}\displaystyle \sum \text{Res}_{A}(G)=-\frac{m_{j}}{|k_{\Vert }|}\unicode[STIX]{x1D6FA}_{j}U\unicode[STIX]{x1D64F}_{n}|_{p_{\Vert }=p_{\text{pole}}}, & & \displaystyle\end{eqnarray}$$

where $p_{\text{pole}}=m_{j}(\unicode[STIX]{x1D714}-n\unicode[STIX]{x1D6FA}_{j})/k_{\Vert }$ is the parallel momentum associated with pole $A$ .

It is a common approach to decompose the background distribution functions in terms of analytical expressions and then to evaluate these at the complex poles. Complete orthogonal basis functions such as Hermite, Legendre or Chebyshev polynomials are the prime candidates for such a decomposition since they can represent $f_{0j}$ to an arbitrary degree of accuracy (Robinson Reference Robinson1990; Weideman Reference Weideman1995; Xie Reference Xie2013). These approaches are useful when $f_{0j}$ deviates only slightly from a Maxwellian. They require, however, very high orders of decomposition and are thus slow in the presence of typical structures that we see in the solar wind such as a proton core–beam configuration. Therefore, they are unsuitable for ALPS’s purpose, and we pursue a different approach, which we call the hybrid analytic continuation. The basic idea behind this approach is to integrate $I$ numerically along the real axis whenever possible and to resort to an algebraic function for the sole purpose of the evaluation of $\text{Res}_{A}(G)$ when necessary.

For the determination of an appropriate algebraic function, ALPS allows the user to choose an arbitrary combination of fit functions to represent $f_{0j}$ and automatically evaluates the fits before the integration begins. The code evaluates the fits separately at each $p_{\bot }$ , so that no assumption is made as to the structure of $f_{0j}$ in the $p_{\bot }$ -direction. ALPS uses these functions only if a pole is within the integration domain and only if $\unicode[STIX]{x1D6FE}\leqslant 0$ . The intrinsic fit functions that the code can combine include a Maxwellian distribution,

(3.10) $$\begin{eqnarray}\displaystyle f_{0j}=\frac{1}{\unicode[STIX]{x03C0}^{3/2}m_{j}^{3}w_{\bot j}^{2}w_{\Vert }}\exp \left(-\frac{p_{\bot }^{2}}{m_{j}^{2}w_{\bot j}^{2}}-\frac{(p_{\Vert }-m_{j}U_{j})^{2}}{m_{j}^{2}w_{\Vert j}^{2}}\right), & & \displaystyle\end{eqnarray}$$

where $w_{\bot j}\equiv \sqrt{2k_{\text{B}}T_{\bot j}/m_{j}}$ ( $w_{\Vert j}\equiv \sqrt{2k_{\text{B}}T_{\Vert j}/m_{j}}$ ) is the thermal speed of species $j$ in the direction perpendicular (parallel) with respect to $\boldsymbol{B}_{0}$ , $T_{\bot j}$ ( $T_{\Vert j}$ ) is the temperature of species $j$ perpendicular (parallel) to $\boldsymbol{B}_{0}$ , $k_{\text{B}}$ is the Boltzmann constant and $U_{j}$ is the $\boldsymbol{B}_{0}$ -parallel drift speed of species $j$ ; a $\unicode[STIX]{x1D705}$ -distribution (Summers et al. Reference Summers, Xue and Thorne1994; Astfalk, Görler & Jenko Reference Astfalk, Görler and Jenko2015),

(3.11) $$\begin{eqnarray}\displaystyle f_{0j} & = & \displaystyle \frac{1}{m_{j}^{3}w_{\bot j}^{2}w_{\Vert j}}\left[\frac{2}{\unicode[STIX]{x03C0}(2\unicode[STIX]{x1D705}-3)}\right]^{3/2}\frac{\tilde{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D705}+1)}{\tilde{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D705}-1/2)}\nonumber\\ \displaystyle & & \displaystyle \times \left\{1+\frac{2}{2\unicode[STIX]{x1D705}-3}\left[\frac{p_{\bot }^{2}}{m_{j}^{2}w_{\bot j}^{2}}+\frac{(p_{\Vert }-m_{j}U_{j})^{2}}{m_{j}^{2}w_{\Vert j}^{2}}\right]\right\}^{-(\unicode[STIX]{x1D705}+1)};\end{eqnarray}$$

and a Jüttner distribution (Jüttner Reference Jüttner1911; Chacón-Acosta, Dagdug & Morales-Técotl Reference Chacón-Acosta, Dagdug and Morales-Técotl2010),

(3.12) $$\begin{eqnarray}\displaystyle f_{0j}=\frac{1}{2\unicode[STIX]{x03C0}m_{j}^{3}cw_{j}^{2}\text{K}_{2}(2c^{2}/w_{j}^{2})}\exp \left(-2\frac{c^{2}}{w_{j}^{2}}\sqrt{1+\frac{|\boldsymbol{p}|^{2}}{m_{j}^{2}c^{2}}}\right); & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is the $\unicode[STIX]{x1D705}$ -index and $\text{K}_{2}$ is the modified Bessel function of the second kind. The Jüttner distribution is the thermodynamic equilibrium distribution if $k_{\text{B}}T_{j}\gtrsim m_{j}c^{2}$ . The exponential in (3.12) reduces to the Maxwellian $\exp (-v^{2}/w_{j}^{2})$ with a different $\boldsymbol{p}$ -independent normalisation factor for $p^{2}/m_{j}^{2}c^{2}\ll 1$ . We use an automated Levenberg–Marquardt-fit algorithm (Levenberg Reference Levenberg1944; Marquardt Reference Marquardt1963) and describe the details of the fit routine in appendix B.

3.3 The poles in a relativistic plasma

The analytic continuation and pole handling in the relativistic case entail a further complication due to the non-trivial $\boldsymbol{p}$ -dependence of the resonant denominator in (2.9) (Buti Reference Buti1962; Lerche Reference Lerche1968). We define a plasma to be relativistic when there is a significant number of particles at relativistic velocities. This can be the case in plasmas with relativistic temperatures ( $k_{\text{B}}T_{j}\gtrsim m_{j}c^{2}$ ) or in plasmas with relativistic beams ( $P_{j}\gtrsim m_{j}c$ , where $P_{j}$ is the drift momentum). Using the relativistic expression for $\unicode[STIX]{x1D6FA}_{j}$ in (2.4) shows that we can write for the pole of the function under the integral sign in (3.1)

(3.13) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D714}-k_{\Vert }v_{\Vert }-n\unicode[STIX]{x1D6FA}_{j}}=-\frac{1}{k_{\Vert }}\frac{\unicode[STIX]{x1D6E4}m_{j}}{\displaystyle \left(p_{\Vert }-\frac{\unicode[STIX]{x1D714}}{k_{\Vert }}\unicode[STIX]{x1D6E4}m_{j}+n\frac{\unicode[STIX]{x1D6FA}_{0j}}{k_{\Vert }}m_{j}\right)}, & & \displaystyle\end{eqnarray}$$

where

(3.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}\equiv \sqrt{1+\frac{p_{\bot }^{2}+p_{\Vert }^{2}}{m_{j}^{2}c^{2}}} & & \displaystyle\end{eqnarray}$$

is the Lorentz factor. We define the dimensionless parallel momentum $\bar{p}_{\Vert }\equiv p_{\Vert }/m_{j}c$ . The dimensionless parallel momentum associated with the relativistic pole is given by

(3.15) $$\begin{eqnarray}\displaystyle \bar{p}_{\text{pole}}=\unicode[STIX]{x1D6E4}\frac{\unicode[STIX]{x1D714}}{k_{\Vert }c}-\frac{n\unicode[STIX]{x1D6FA}_{0j}}{k_{\Vert }c}. & & \displaystyle\end{eqnarray}$$

We apply the technique proposed by Lerche (Reference Lerche1967) to transform (2.9) from the $(p_{\bot },p_{\Vert })$ coordinate system to the $(\unicode[STIX]{x1D6E4},\bar{p}_{\Vert })$ coordinate system (see also Swanson Reference Swanson2002; Lazar & Schlickeiser Reference Lazar and Schlickeiser2006; López et al. Reference López, Moya, Muñoz, Viñas and Valdivia2014, Reference López, Moya, Navarro, Araneda, Muñoz, Viñas and Valdivia2016). This transformation yields

(3.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D74C}_{j} & = & \displaystyle 2\unicode[STIX]{x03C0}m_{j}^{3}c^{3}\frac{\unicode[STIX]{x1D714}_{\text{p}j}^{2}}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FA}_{0j}}\int _{1}^{\infty }\text{d}\unicode[STIX]{x1D6E4}\int _{-\sqrt{\unicode[STIX]{x1D6E4}^{2}-1}}^{+\sqrt{\unicode[STIX]{x1D6E4}^{2}-1}}\text{d}\bar{p}_{\Vert }\left[\hat{\boldsymbol{e}}_{\Vert }\hat{\boldsymbol{e}}_{\Vert }\frac{\unicode[STIX]{x1D6FA}_{0j}}{\unicode[STIX]{x1D714}}\bar{p}_{\Vert }\frac{\unicode[STIX]{x2202}f_{0j}}{\unicode[STIX]{x2202}\bar{p}_{\Vert }}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,\mathop{\sum }_{n=-\infty }^{+\infty }\frac{\unicode[STIX]{x1D6FA}_{0j}}{k_{\Vert }c}\left(\frac{\unicode[STIX]{x2202}f_{0j}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}+\frac{k_{\Vert }c}{\unicode[STIX]{x1D714}}\frac{\unicode[STIX]{x2202}f_{0j}}{\unicode[STIX]{x2202}\bar{p}_{\Vert }}\right)\frac{1}{\displaystyle \bar{p}_{\Vert }-\unicode[STIX]{x1D6E4}\frac{\unicode[STIX]{x1D714}}{k_{\Vert }c}+\frac{n\unicode[STIX]{x1D6FA}_{0j}}{k_{\Vert }c}}\bar{\unicode[STIX]{x1D64F}}_{n}\right],\end{eqnarray}$$

where

(3.17) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D64F}_{n}}\equiv \left(\begin{array}{@{}ccc@{}}\displaystyle \frac{n^{2}\text{J}_{n}^{2}}{\bar{z}^{2}} & \displaystyle \frac{\text{i}n\text{J}_{n}\text{J}_{n}^{\prime }}{\bar{z}}\bar{p}_{\bot } & \displaystyle \frac{n\text{J}_{n}^{2}\bar{p}_{\Vert }}{\bar{z}}\\ \displaystyle -\frac{\text{i}n\text{J}_{n}\text{J}_{n}^{\prime }}{\bar{z}}\bar{p}_{\bot } & \displaystyle (\text{J}_{n}^{\prime })^{2}\bar{p}_{\bot }^{2} & \displaystyle -\text{i}\text{J}_{n}\text{J}_{n}^{\prime }\bar{p}_{\Vert }\bar{p}_{\bot }\\ \displaystyle \frac{n\text{J}_{n}^{2}\bar{p}_{\Vert }}{\bar{z}} & \displaystyle \text{i}\text{J}_{n}\text{J}_{n}^{\prime }\bar{p}_{\Vert }\bar{p}_{\bot } & \displaystyle \text{J}_{n}^{2}\bar{p}_{\Vert }^{2}\end{array}\right), & & \displaystyle\end{eqnarray}$$

$\bar{p}_{\bot }\equiv \sqrt{\unicode[STIX]{x1D6E4}^{2}-1-\bar{p}_{\Vert }^{2}}$ , $\bar{z}\equiv k_{\bot }c/\unicode[STIX]{x1D6FA}_{0j}$ and the Bessel functions are evaluated as $\text{J}_{n}\equiv \text{J}_{n}(\bar{z}\bar{p}_{\bot })$ . Whenever ALPS performs a relativistic calculation and

(3.18) $$\begin{eqnarray}\displaystyle -P_{\max ,\Vert j}\leqslant \text{Re}(\bar{p}_{\text{pole}})\leqslant +P_{\max ,\Vert j}, & & \displaystyle\end{eqnarray}$$

the code automatically transforms from $(p_{\bot },p_{\Vert })$ to $(\unicode[STIX]{x1D6E4},\bar{p}_{\Vert })$ coordinates and applies the polyharmonic-spline algorithm described in appendix C to create an equally spaced and homogeneous grid in $(\unicode[STIX]{x1D6E4},\bar{p}_{\Vert })$ coordinates. In this coordinate system, we perform the integration near poles and the analytic continuation in the same way as described in §§ 3.1 and 3.2, but using the relativistic parallel momentum associated with the pole from (3.15). For reasons of numerical performance, we use the integration based on (3.16) only if there is a pole within the integration domain. Otherwise, we employ the faster integration method based on (2.9) even in the relativistic case.

4.1 Maxwellian distributions

There are numerous codes for the hot-plasma dispersion relation in a plasma with Maxwellian or bi-Maxwellian background distributions. We use our code PLUME (Klein & Howes Reference Klein and Howes2015) for an electron–proton plasma and calculate the dispersion relations of Alfvén/ion-cyclotron (A/IC) and fast-magnetosonic/whistler (FM/W) waves. We then set up Maxwellian $f_{0}$ tables with the same parameters as those used with PLUME and calculate the dispersion relations based on these $f_{0}$ tables with ALPS. We compare the PLUME and ALPS results for quasi-parallel and quasi-perpendicular propagation in figure 1. The panels show both the real part of the frequency $\unicode[STIX]{x1D714}_{\text{r}}$ and its imaginary part $\unicode[STIX]{x1D6FE}$ as functions of the parallel and perpendicular wavenumbers, respectively. We use $\unicode[STIX]{x1D6FD}_{\bot j}=\unicode[STIX]{x1D6FD}_{\Vert j}=1$ for both protons and electrons, and $v_{\text{A}}/c=10^{-4}$ , where $\unicode[STIX]{x1D6FD}_{\bot j}\equiv 8\unicode[STIX]{x03C0}n_{0j}k_{\text{B}}T_{\bot j}/B_{0}^{2}$ and $\unicode[STIX]{x1D6FD}_{\Vert j}\equiv 8\unicode[STIX]{x03C0}n_{0j}k_{\text{B}}T_{\Vert j}/B_{0}^{2}$ . We normalise all frequencies in units of the proton-cyclotron frequency $\unicode[STIX]{x1D6FA}_{0\text{p}}$ and all length scales in units of the proton skin depth $d_{\text{p}}\equiv v_{\text{A}}/\unicode[STIX]{x1D6FA}_{0\text{p}}$ . The momentum-space resolution for the ALPS calculation in the quasi-parallel limit is $N_{\bot }=320$ , $N_{\Vert }=640$ , $P_{\max ,\Vert \text{p}}=8m_{\text{p}}v_{\text{A}}$ and $P_{\max ,\Vert \text{e}}=0.19m_{\text{p}}v_{\text{A}}$ . In the quasi-perpendicular limit, we use $N_{\bot }=240$ , $N_{\Vert }=480$ , $P_{\max ,\Vert \text{p}}=6m_{\text{p}}v_{\text{A}}$ and $P_{\max ,\Vert \text{e}}=0.14m_{\text{p}}v_{\text{A}}$ . In both cases, we set $P_{\max ,\bot j}=P_{\max ,\Vert j}$ , $J_{\max }=10^{-45}$ , $M_{\text{I}}=5$ , $M_{\text{P}}=100$ , $t_{\text{lim}}=0.01$ and scan through $\boldsymbol{k}$ -space in 256 logarithmically spaced steps from $k_{\Vert }d_{\text{p}}=10^{-2}$ to $k_{\Vert }d_{\text{p}}=10$ and $k_{\bot }d_{\text{p}}=10^{-2}$ to $k_{\bot }d_{\text{p}}=10$ , respectively. We study the accuracy of the results depending on the resolution in § A.1. Figure 1 shows that ALPS reproduces these Maxwellian examples very well. We note that these plasma parameters represent typical solar-wind conditions at 1 au.

Figure 1.

Dispersion relations for the A/IC wave and the FM/W wave in a Maxwellian plasma in quasi-parallel (a) and quasi-perpendicular (b) propagation. For the calculations shown on (a), we keep $k_{\bot }d_{\text{p}}=10^{-3}$ constant and scan through $k_{\Vert }$ . For the calculations shown on (b), we keep $k_{\Vert }d_{\text{p}}=10^{-3}$ constant and scan through $k_{\bot }$ . The A/IC mode in quasi-perpendicular propagation corresponds to the kinetic Alfvén wave (KAW) at $k_{\bot }d_{\text{p}}\gtrsim 1/\sqrt{\unicode[STIX]{x1D6FD}_{\Vert \text{p}}}$ . We compare ALPS with the standard Maxwellian solutions from PLUME for an electron–proton plasma with the same plasma parameters. Both numerical models agree well in both the real part $\unicode[STIX]{x1D714}_{\text{r}}$ of the frequency and its imaginary part $\unicode[STIX]{x1D6FE}$ .

Figure 2.

Comparison of dispersion maps from PLUME (a) and ALPS (b) for $k_{\bot }d_{\text{p}}=k_{\Vert }d_{\text{p}}=10^{-3}$ . The lines show isocontours of constant $\lg |\text{det}\,{\mathcal{D}}|$ . The colour scheme varies from small $\lg |\text{det}\,{\mathcal{D}}|$ in blue to large $\lg |\text{det}\,{\mathcal{D}}|$ in red. Minima in this map correspond to solutions to the hot-plasma dispersion relation.

In order to illustrate another representation of the plasma dispersion relation, we show a comparison of dispersion maps from PLUME and ALPS in figure 2. Dispersion maps are diagrams of isocontours of constant $\lg |\text{det}\,{\mathcal{D}}|$ , where ${\mathcal{D}}$ is the tensor from (2.12), in the $\unicode[STIX]{x1D714}_{\text{r}}$ – $\unicode[STIX]{x1D6FE}$ plane. Solutions to the hot-plasma dispersion relation appear as local minima in these diagrams. After calculating a dispersion map for a given $\boldsymbol{k}$ , ALPS automatically determines the locations of its local minima in terms of $\unicode[STIX]{x1D714}_{\text{r}}$ and $\unicode[STIX]{x1D6FE}$ . These minima represent ideal initial guesses for the Newton-secant root-finding search when scanning through $\boldsymbol{k}$ . Although the calculation of a dispersion map still requires the calculation of all $\unicode[STIX]{x1D74C}_{j}$ , it does not entail the application of the Newton-secant root-finding algorithm itself. We use a Maxwellian plasma model with $\unicode[STIX]{x1D6FD}_{\bot j}=\unicode[STIX]{x1D6FD}_{\Vert j}=1$ for both protons and electrons, $k_{\bot }d_{\text{p}}=k_{\Vert }d_{\text{p}}=10^{-3}$ and $v_{\text{A}}/c=10^{-4}$ . For the ALPS calculation, we use $N_{\bot }=240$ , $N_{\Vert }=480$ , $P_{\max ,\Vert \text{p}}=P_{\max ,\bot \text{p}}=6m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\Vert \text{e}}=P_{\max ,\bot \text{e}}=0.14m_{\text{p}}v_{\text{A}}$ , $J_{\max }=10^{-45}$ , $M_{\text{I}}=5$ , $M_{\text{P}}=100$ and $t_{\text{lim}}=0.01$ . Both the PLUME and the ALPS calculations reveal seven solutions to the dispersion relation. We note that the point $\unicode[STIX]{x1D714}_{\text{r}}=\unicode[STIX]{x1D6FE}=0$ is a maximum and does not represent a solution to the dispersion relation. The solutions at $\unicode[STIX]{x1D714}_{\text{r}}=\pm 10^{-3}\unicode[STIX]{x1D6FA}_{0\text{p}}$ and $\unicode[STIX]{x1D6FE}=-2.3\times 10^{-10}\unicode[STIX]{x1D6FA}_{0\text{p}}$ are the forward and backward propagating A/IC waves. The solutions at $\unicode[STIX]{x1D714}_{\text{r}}=\pm 2\times 10^{-3}\unicode[STIX]{x1D6FA}_{0\text{p}}$ and $\unicode[STIX]{x1D6FE}=-5.4\times 10^{-5}\unicode[STIX]{x1D6FA}_{0\text{p}}$ are the forward and backward propagating FM/W waves. The solutions at $\unicode[STIX]{x1D714}_{\text{r}}=\pm 1.2\times 10^{-3}\unicode[STIX]{x1D6FA}_{0\text{p}}$ and $\unicode[STIX]{x1D6FE}=-7.3\times 10^{-4}\unicode[STIX]{x1D6FA}_{0\text{p}}$ are the forward and backward propagating slow waves (ion-acoustic waves). Lastly, the solution at $\unicode[STIX]{x1D714}_{\text{r}}=0$ and $\unicode[STIX]{x1D6FE}=-7.2\times 10^{-4}\unicode[STIX]{x1D6FA}_{0\text{p}}$ is the non-propagating slow mode, which is sometimes denoted ‘entropy mode’, (Verscharen et al. Reference Verscharen, Chandran, Klein and Quataert2016; Verscharen, Chen & Wicks Reference Verscharen, Chen and Wicks2017). The comparison of both panels in figure 2 shows that ALPS reproduces these seven plasma modes under typical solar-wind conditions in the Maxwellian limit.

Figure 3.

Comparison of dispersion relations for the A/IC instability (a) and the mirror-mode instability (b) from PLUME and ALPS. We use $T_{\bot \text{p}}/T_{\Vert \text{p}}=3$ . For the calculation of the A/IC instability, we keep $k_{\bot }d_{\text{p}}=10^{-3}$ constant and scan through $k_{\Vert }$ . For the calculation of the mirror-mode instability, we keep the angle $\unicode[STIX]{x1D703}=75^{\circ }$ constant and scan through $|\boldsymbol{k}|$ .

4.2 Anisotropic bi-Maxwellian distributions

PLUME, like most other standard hot-plasma dispersion-relation solvers, allows us to use anisotropic bi-Maxwellian representations for the background distribution functions. Such a configuration can lead to instability if the temperature anisotropy exceeds the threshold for an anisotropy-driven plasma instability. As an example for a propagating instability, we calculate the dispersion relation for the parallel A/IC instability (Harris Reference Harris1961; Davidson & Ogden Reference Davidson and Ogden1975; Yoon et al. Reference Yoon, Seough, Khim, Kim, Kwon, Park, Parkh and Park2010), and as an example for a non-propagating instability, we calculate the dispersion relation for the mirror-mode instability (Rudakov & Sagdeev Reference Rudakov and Sagdeev1961; Tajiri Reference Tajiri1967; Southwood & Kivelson Reference Southwood and Kivelson1993). The thresholds for both of these instabilities fulfil $T_{\bot \text{p}}>T_{\Vert \text{p}}$ . For this demonstration, we use PLUME to calculate $\unicode[STIX]{x1D714}_{\text{r}}$ and $\unicode[STIX]{x1D6FE}$ as functions of the wavenumber in a plasma with bi-Maxwellian protons and Maxwellian electrons using $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=\unicode[STIX]{x1D6FD}_{\bot \text{e}}=1$ , $T_{\bot \text{p}}/T_{\Vert \text{p}}=3$ and $v_{\text{A}}/c=10^{-4}$ . We then set-up bi-Maxwellian $f_{0}$ tables with the same parameters and calculate the dispersion relations for both instabilities with ALPS. We show the results in figure 3. For the ALPS calculation, we use $N_{\bot }=320$ , $N_{\Vert }=640$ , $P_{\max ,\Vert \text{p}}=8m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\bot \text{p}}=13.9m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\Vert \text{e}}=P_{\max ,\bot \text{e}}=0.19m_{\text{p}}v_{\text{A}}$ , $J_{\max }=10^{-45}$ , $M_{\text{I}}=5$ , $M_{\text{P}}=100$ and $t_{\text{lim}}=0.01$ . For the A/IC instability, we scan through $k_{\Vert }$ from $k_{\Vert }d_{\text{p}}=10^{-2}$ to $k_{\Vert }d_{\text{p}}=10$ in 256 logarithmically spaced steps. For the mirror-mode instability, we scan through $k$ from $kd_{\text{p}}=10^{-2}$ to $kd_{\text{p}}=5$ in 154 logarithmically spaced steps and keep the angle $\unicode[STIX]{x1D703}\equiv \text{arcsin}(k_{\bot }/k)=75^{\circ }$ constant. We study the accuracy of these results depending on the resolution in § A.2.

Both PLUME and ALPS show that the A/IC wave and the mirror mode are unstable in different wave-vector ranges for the given parameter set. The good agreement between the PLUME solutions and the ALPS solutions shows that ALPS successfully calculates the dispersion relations of both instabilities in a bi-Maxwellian plasma. Although the correct mirror mode exhibits $\unicode[STIX]{x1D714}_{\text{r}}=0$ , the ALPS solution shows a very small yet finite $\unicode[STIX]{x1D714}_{\text{r}}$ . We discuss this error and its dependence on the velocity-space resolution in § A.2.

4.3 Anisotropic $\unicode[STIX]{x1D705}$ -distributions

Astfalk et al. (Reference Astfalk, Görler and Jenko2015) developed the code DSHARK to calculate dispersion relations in plasmas with bi- $\unicode[STIX]{x1D705}$ -distributions. As one example, these authors discuss the FM/W instability in an anisotropic electron–proton plasma with $\unicode[STIX]{x1D705}_{\text{p}}=\unicode[STIX]{x1D705}_{\text{e}}=8$ , $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=2$ , $\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=4$ , $T_{\bot \text{p}}/T_{\Vert \text{p}}=0.4$ and $T_{\bot \text{e}}/T_{\Vert \text{e}}=0.5$ (see figure 1 from Astfalk et al. Reference Astfalk, Görler and Jenko2015). The angle between $\boldsymbol{k}$ and $\boldsymbol{B}_{0}$ is constant for this calculation and set to $\unicode[STIX]{x1D703}=0.001^{\circ }$ . We use DSHARK to reproduce this test case and set up $\unicode[STIX]{x1D705}$ -distributed $f_{0}$ tables with the same parameters in order to compare the DSHARK results with ALPS. We show this comparison in figure 4. In ALPS, we use $N_{\bot }=400$ , $N_{\Vert }=800$ , $P_{\max ,\Vert \text{p}}=10m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\bot \text{p}}=6.32m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\Vert \text{e}}=0.33m_{\text{p}}v_{\text{A}}$ , $P_{\max ,\bot \text{e}}=0.23m_{\text{p}}v_{\text{A}}$ , $J_{\max }=10^{-45}$ , $M_{\text{I}}=5$ , $M_{\text{P}}=500$ , $t_{\text{lim}}=0.01$ and $v_{\text{A}}/c=10^{-4}$ . We scan through $k_{\Vert }$ from $k_{\Vert }d_{p\text{}}=10^{-1}$ to $k_{\Vert }d_{\text{p}}\approx 3.42$ , where $\unicode[STIX]{x1D714}_{\text{r}}$ changes its sign after 87 logarithmically spaced steps.

Figure 4.

Comparison of dispersion relations of the FM/W instability in a $\unicode[STIX]{x1D705}$ -distributed plasma from DSHARK and ALPS. In both plasma models, we keep $\unicode[STIX]{x1D703}=0.001^{\circ }$ constant and scan through $k_{\Vert }$ . (a) shows the real part of the wave frequency, and (b) shows its imaginary part.

ALPS reproduces the DSHARK results for the FM/W instability well. The results also agree with the previous work by Lazar et al. (Reference Lazar, Poedts and Schlickeiser2011).

4.4 Relativistic Jüttner distributions

As one example for a dispersion relation in a relativistic plasma, we reproduce the results by López et al. (Reference López, Moya, Muñoz, Viñas and Valdivia2014) for an electron–positron pair plasma with a Jüttner distribution using $v_{\text{A}}/c=1$ , $m_{\text{p}}=m_{\text{e}}$ , $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=(0.2,0.4,1.0)$ and $T_{\bot j}=T_{\Vert j}$ for both positrons and electrons. We set up a Jüttner-distributed $f_{0}$ table with the same parameters and calculate the dispersion relations of the A/IC wave and the Ordinary wave (O-mode) in the plasma, keeping the perpendicular wavenumber constant at $k_{\bot }d_{\text{p}}=10^{-3}$ . We use $N_{\bot }=30$ , $N_{\Vert }=60$ , $P_{\max }=5m_{\text{p}}v_{\text{A}}$ , $J_{\max }=10^{-45}$ , $M_{\text{I}}=5$ , $M_{\text{P}}=300$ and $t_{\text{lim}}=0.01$ . Our interpolation method transforms the $(p_{\bot },p_{\Vert })$ grid to the $(\unicode[STIX]{x1D6E4},\bar{p}_{\Vert })$ grid with $N_{\unicode[STIX]{x1D6E4}}=500$ and $N_{\bar{p}_{\Vert }}=500$ steps in $\unicode[STIX]{x1D6E4}$ and $\bar{p}_{\Vert }$ , respectively. We scan through $k_{\Vert }$ from $k_{\Vert }d_{\text{p}}=10^{-1}$ to $k_{\Vert }d_{\text{p}}=10$ in 128 logarithmically spaced steps. We show the results in figure 5. López et al. (Reference López, Moya, Muñoz, Viñas and Valdivia2014) show their results for these parameters in their figure 1. Our comparison with the ALPS dispersion relation in figure 5 shows a good agreement and confirms our relativistic model. The deviation between the results from López et al. (Reference López, Moya, Muñoz, Viñas and Valdivia2014) and ALPS is only visible in the real part of the frequency at the large- $k_{\Vert }$ /low- $\unicode[STIX]{x1D714}_{\text{r}}$ end of the A/IC branches.

Figure 5.

Dispersion relations of the quasi-parallel A/IC wave (solutions at low $\unicode[STIX]{x1D714}_{\text{r}}$ ) and the ordinary wave (solutions at high $\unicode[STIX]{x1D714}_{\text{r}}$ ) in a relativistic electron–positron pair plasma with Jüttner distributions. We keep $k_{\bot }d_{\text{p}}=10^{-3}$ constant and scan through $k_{\Vert }$ . (a) shows the real part of the frequency, and (b) shows the imaginary part of the frequency. The lines show ALPS solutions, and the crosses show the results from figure 1 of López et al. (Reference López, Moya, Muñoz, Viñas and Valdivia2014). The three colours correspond to $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=0.2$ (red), $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=0.4$ (green) and $\unicode[STIX]{x1D6FD}_{\Vert \text{p}}=\unicode[STIX]{x1D6FD}_{\Vert \text{e}}=1.0$ (blue) in both panels and for both modes.