ALPS: The Arbitrary Linear Plasma Solver

The Arbitrary Linear Plasma Solver (ALPS) is a parallelised numerical code that solves the dispersion relation in a hot (even relativistic) magnetised plasma with an arbitrary number of particle species with arbitrary gyrotropic equilibrium distribution functions for any direction of wave propagation with respect to the background field. ALPS reads the background momentum distributions as tables of values on a $(p_{\perp},p_{\parallel})$ grid, where $p_{\perp}$ and $p_{\parallel }$ are the momentum coordinates in the directions perpendicular and parallel to the background magnetic field, respectively. We present the mathematical and numerical approach used by ALPS and introduce our algorithms for the handling of poles and the analytic continuation for the Landau contour integral. We then show test calculations of dispersion relations for a selection of stable and unstable configurations in Maxwellian, bi-Maxwellian, $\kappa$-distributed, and J\"uttner-distributed plasmas. These tests demonstrate that ALPS derives reliable plasma dispersion relations. ALPS will make it possible to determine the properties of waves and instabilities in the non-equilibrium plasmas that are frequently found in space, laboratory experiments, and numerical simulations.


Introduction
The vast majority of the visible matter in the universe is in the plasma state. The solar wind is an example of such an astrophysical plasma. Due to its accessibility to spacecraft, it is the perfect environment for making comparisons between theoretical plasma-physics predictions and in-situ observations in the astrophysical context with access to wide scale separations (see, for example, Marsch 2006). Plasma can deviate from thermodynamic equilibrium if the relaxation due to particle collisions occurs on timescales that are larger than the characteristic timescales of the collective plasma behaviour. Such a collisionless plasma is characterised by non-Maxwellian features in its velocity distribution functions. In the fast solar wind, this condition is frequently fulfilled, and, consequently, the observed distribution functions often deviate from the entropically favoured Maxwellian shape (Vasyliunas 1968;Gosling et al. 1981;Lui & Krimigis 1981;Marsch et al. 1982b,a;Armstrong et al. 1983;Lui & Krimigis 1983;Christon et al. 1988;Williams et al. 1988). In particular, beams and temperature anisotropies are some of the observed features in the distributions of ions and electrons in the solar wind (Pilipp et al. 1987a,b;Hellinger et al. 2006;Marsch 2006;Bale et al. 2009). If these deviations from equilibrium are suitably extreme, the plasma becomes unstable and generates waves or non-propagating structures that react back upon the plasma to reduce the deviations from equilibrium (Eviatar & Schulz 1970;Schwartz 1980;Gary 1993;Hellinger & Trávníček 2011, 2013. The behaviour of plasma waves and instabilities is typically studied with the help of numerical codes that solve the hot-plasma dispersion relation. Traditionally, these codes (like WHAMP, PLUME, or NHDS) use a shifted bi-Maxwellian background distribution function as the zeroth-order description for the plasma state (Roennmark 1982;Quataert 1998;Klein et al. 2012;Verscharen et al. 2013a). For nearly collisionless plasmas, however, the bi-Maxwellian distribution function is a mathematical convenience rather than a reliable representation of the true plasma distribution function, and many space-plasma observations show that the bi-Maxwellian representation is not accurate (Hundhausen 1970;Leubner 1978;Marsch et al. 1982b;Pilipp et al. 1987a;Marsch & Tu 2001;Štverák et al. 2009). Some previous approaches in non-Maxwellian solvers treated certain limits or geometries (Dum et al. 1980;Summers & Thorne 1991;Summers et al. 1994;Xue et al. 1993Xue et al. , 1996Hellberg et al. 2005;Cattaert et al. 2007;Lazar, M. & Poedts, S. 2009;Mace & Sydora 2010;Lazar et al. 2011;Galvaõ et al. 2012;Xie 2013;Lazar & Poedts 2014; or faced challenges in the weakly-damped limit (Hellinger & Trávníček 2011).
We present our numerical code ALPS (Arbitrary Linear Plasma Solver), which solves the full hot-plasma dispersion relation in a plasma consisting of an arbitrary number of particle species with arbitrary background distribution functions f 0j and with arbitrary directions of wave propagation with respect to the uniform background magnetic field. ALPS is also able to solve the dispersion relation for relativistic plasmas. Matsuda & Smith (1992) developed a code similar to ALPS that calculates the dispersion relation in an arbitrary plasma with relativistic effects. Their code uses a cubic spline fit to both fill data gaps and approximate the analytic continuation, while ALPS uses a novel method called hybrid analytic continuation. The spline method forfeits its accuracy for strongly damped solutions since the calculation of the dispersion relation requires the evaluation of the spline at a complex value that is distant from the real grid points by which the spline is supported. Our method does not suffer from this problem. Astfalk & Jenko (2017) also use a cubic-spline interpolation for the analytic continuation and as the basis for the integration in their code LEOPARD. This procedure allows for algebraic simplifications that enhance the speed of the integration significantly. LEOPARD, however, does not capture relativistic effects.
In Section 2, we review the underlying theory of the hot-plasma dispersion relation. Section 3 presents ALPS's numerical approach. In Section 4, we compare ALPS results to known limits of the hot-plasma dispersion relation such as Maxwellian, bi-Maxwellian, κ-distributed, and relativistic pair plasmas. In Section 5, we discuss our results and the applicability of ALPS to measured plasma distributions. The Appendix describes how ALPS solutions depend on the resolution of the background distributions, discusses of the Levenberg-Marquardt-fit routine used in our hybrid-analytic-continuation method, and describes our strategy for numerically refining coarse-grained distribution functions obtained from spacecraft measurements.

The Linear Dispersion Relation of a Hot Plasma
In this section, we discuss the mathematical basis for the calculation of the hotplasma dispersion relation following the presentation and notation of Stix (1992). The determination of the kinetic wave dispersion relation in a hot plasma is based on the linearised set of Maxwell's equations and the linearised Vlasov equation (Stix 1992;Gary 1993). A wave or instability is then associated with a first-order perturbation δf j in the distribution function of species j about a prescribed time-averaged background distribution function f 0j , where x is the spatial coordinate and p is the momentum coordinate. As with the distribution function f j in Equation (2.1), we take the magnetic field B to be the sum of a uniform background magnetic field B 0 and a fluctuating magnetic field δB. We assume that E = δE; i.e., the average electric field is zero. Linear theory expresses δf j as a function of f 0j and the electromagnetic field components.
The distribution function f j in a collisionless plasma evolves according to the Vlasov equation, where q j is the charge of a particle of species j, c is the speed of light, and v is the velocity coordinate. We assume that all fluctuating quantities behave like plane waves; i.e., ∝ exp (ik · x − iωt), where k is the wave vector and ω is the (complex) frequency. Linearising Equation (2.2), using Faraday's law, and applying the method of characteristics, we obtain where E = (E x , E y , E z ) is the electric field, φ is the azimuthal angle of the momentum vector p, ϑ is the azimuthal angle of the wavevector k, and the index ⊥ ( ) refers to the direction perpendicular (parallel) with respect to the background magnetic field B 0 , is the relativistic gyrofrequency, m j is the rest mass of a particle of species j, (2.7) The first velocity moments of the distribution functions of all species define the current density j through where χ j is the contribution of species j to the plasma susceptibility. Without loss of generality, we choose a cylindrical coordinate system in which k y = ϑ = 0 and apply a set of Bessel-function identities in order to facilitate the integration over φ and τ in Equation (2.3). This allows us to rewrite the plasma susceptibilities as (provided that Im(ω) > 0) where ω pj ≡ 4πn j q 2 j /m j is the plasma frequency of species j, Ω 0j ≡ q j B 0 /m j c is the non-relativistic gyrofrequency, n j is the density of species j, and the tensor T n is defined as where z ≡ k ⊥ v ⊥ /Ω j , and J n ≡ J n (z) is the nth-order Bessel function. For Im(ω) 0, the integral over p is executed as the Landau integral after analytic continuation (for details, see Chapt. 8 of Stix 1992). Equation (2.9) describes the susceptibility for a general background distribution function f 0j in a relativistic plasma. The only assumptions are gyrotropy in f 0j and small amplitudes in the fluctuations so that linearisation is applicable, and a uniform, stationary equilibrium. The numerical challenge in the solution of the plasma dispersion relation results from the integrals over p ⊥ and p in Equation (2.9). We note that, in numerous classical codes for calculation of the linear hot-plasma dispersion relation (Roennmark 1982;Gary 1993;Verscharen et al. 2013b;Klein & Howes 2015), these integrals are greatly simplified by assuming that f 0j is a (bi)-Maxwellian. The dielectric tensor ε of the plasma is related to the plasma susceptibilities from Equation (2.9) through (2.11) Finally, combining Faraday's law and Ampère's law leads to the wave equation, where n ≡ kc/ω is the index . By setting det D = 0, we obtain the dispersion relations ω = ω(k) for non-trivial solutions to Equation (2.12). We write these solutions in the form ω = ω r + iγ, where ω r = Re(ω) and γ = Im(ω).

Numerical Approach
In order to find the solutions to the hot-plasma dispersion relation, ALPS determines the values of ω r and γ that solve Equation (2.12) for specified background distributions f 0j at a given set of values for k, m j , q j , n j , and v A /c, where v A ≡ B 0 / 4πn p m p . ALPS uses an efficient iterative Newton-secant algorithm to solve Equation (2.12) based on an initial guess for ω r and γ (Press et al. 1992). The numerically challenging part for this calculation is the evaluation of χ j in Eqs. (2.9). In the following, we present ALPS's strategy for this evaluation in the non-relativistic case. We discuss the extension to relativistic cases with poles in the integration domain in Section 3.3, which is equivalent to the non-relativistic case with the exception that the coordinate system is transformed from (p ⊥ , p ) to (Γ,p ) and that Equation (3.16) below is used instead of Equation (2.9).
We prescribe the shape of f 0j in input files for each species (called "f 0 table") as an ASCII table that lists p ⊥ , p , and the associated values of f 0j . From this table, we calculate ∂f 0j /∂p ⊥ and ∂f 0j /∂p on the same grid as the f 0 table using second-order finite differencing. The resolution of the f 0 table is given by n ⊥ points in the p ⊥ -direction and n points in the p -direction. The table spans from p ⊥ = 0 to p ⊥ = P max,⊥j in the perpendicular direction and from p = −P max, j to p = P max, j in the parallel direction.
The integration in Equation (2.9) allows us to integrate separately and independently for each n and j. This provides us with a very natural way to parallelise the calculation scheme by assigning the separate integrations to different processors. We use MPI for the parallelisation. The integrating nodes return their contributions to χ j to the master node, which then sums up the contributions, determines the value of ε, and updates the values of ω r and γ through a Newton-secant step. The updated values for ω r and γ are then returned to the integrating nodes, which afterwards evaluate the integration of their updated contribution to χ j . We evaluate all values of n up to a value of ±n max , which is determined as the value of n for which the maximum value of |J n | is smaller than the user-defined parameter J max . The necessary value of n max depends on the wavenumber, the direction of propagation of the treated wave, and the thermal speeds of the plasma components. In bi-Maxwellian codes under typical solar-wind conditions, the accuracy of the dispersion relation is better than ∆|ω|/|ω| ∼ 10 −5 for J max ∼ 10 −45 (which typically corresponds to n max 10 at proton scales).
We use a standard two-dimensional trapezoidal integration scheme to integrate over p ⊥ and p . However, this scheme breaks down near the poles of the integrand in Equation (2.9) and requires a special treatment of the analytic continuation when γ 0. In the remainder of this section, we discuss our strategies to resolve these numerical difficulties.

Integrating Near Poles
A challenge concerning the numerical integration is the treatment of the poles that occur in the term proportional to T n in Equation (2.9). The integrals in question are of the form for γ > 0. For sufficiently small γ, the denominator in Equation (3.1) can become very small along the real p 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 Equation (3.1) in the more generic form where x, t r , and t i are real, g(x) is a smooth function, and the integration is performed along the real axis. We choose a symmetric interval [t r − ∆, t r + ∆] around t r where ∆ g(t r )/g (t r ), and write where "rest" refers to the integration outside the . Following Longman (1958) and Davis & Rabinowitz (1984), we then separate the integrand into its odd and even parts with respect to t r as The integrand in the first integral in Equation (3.4) is odd with respect to t r and thus vanishes after the integration over the symmetric interval around t r . The second integral, on the other hand, is even with respect to t r and thus We define ∆ through a user-defined parameter n I so that ∆ ≡ n I ∆p . We then define δ ≡ ∆/n P , where n P is another user-defined parameter. Except for cases in which |t i | is extremely small, we apply a trapezoidal integration over n P steps of width δ to the integral in Equation (3.5). The smoothness of g(x) allows us to expand g(x) around the nearest grid point of the n I grid points in the interval [t r , t r + ∆] using a Taylor series. By taking ∆p to be sufficiently small, we can retain just the first two terms in the series without losing significant accuracy. Since the integral in Equation (3.5) does not converge numerically if |t i | is extremely small, we implement the following procedure when |t i | t lim , where t lim is a user-defined parameter. We first rewrite Equation (3.5) using truncated Taylor expansions of g(x) and g(2t r − x) around x = t r as We determine g(t r ) and g (t r ) through linear interpolation between the neighbouring grid points to t r . The term proportional to g (t r ) in Equation (3.6) converges numerically for any value of t i . We set the term proportional to g(t r ) equal to its small-t i limit, namely We use this method for both the integration of χ j near poles and the principal-value integration that is necessary if γ = 0.

Analytic Continuation
If γ 0, the integration in Equation (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 ⊥ , p ) at a complex value for p in the nonrelativistic case. In our case, however, f 0j is only defined on a real grid in p ⊥ and p , yet the analytic continuation of f 0j is still uniquely defined. This leads to the known mathematical problem of numerical analytic continuation (Cannon & Miller 1965;Reichel 1986;Fujiwara et al. 2007;Fu et al. 2012;Zhang Z.-Q. Ma 2013;Kranich 2014). Our solution for this problem is our hybrid analytic continuation scheme. We note that this approach is only relevant for damped modes, i.e., γ 0.
Landau's rule of integration around singularities (Landau 1946;Lifshitz & Pitaevskii 1981) leads to the following three cases with the appropriate residues for the evaluation of I(p ⊥ ) for general γ: where C 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 Equation (3.8) are performed along the real axis, and 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 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 1990;Weideman 1995;Xie 2013). 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 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 ⊥ , so that no assumption is made as to the structure of f 0j in the p ⊥ -direction. ALPS uses these functions only if a pole is within the integration domain and only if γ 0. The intrinsic fit functions that the code can combine include a Maxwellian distribution, is the thermal speed of species j in the direction perpendicular (parallel) with respect to B 0 , T ⊥j (T j ) is the temperature of species j perpendicular (parallel) to B 0 , k B is the Boltzmann constant, and U j is the B 0 -parallel drift speed of species j; a κ-distribution (Summers et al. 1994;Astfalk et al. 2015), and a Jüttner distribution (Jüttner 1911;Chacón-Acosta et al. 2010), (3.12) where κ is the κ-index,Γ is the gamma function, and K 2 is the modified Bessel function of the second kind. The Jüttner distribution is the thermodynamic-equilibrium distribution if k B T j m j c 2 . The exponential in Equation (3.12) reduces to the Maxwellian exp(−v 2 /w 2 j ) with a different p-independent normalisation factor for p 2 /m 2 j c 2 1. We use an automated Levenberg-Marquardt-fit algorithm (Levenberg 1944;Marquardt 1963) and describe the details of the fit routine in Appendix B.

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 p-dependence of the resonant denominator in Equation (2.9) (Buti 1962;Lerche 1968). 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 B T j m j c 2 ) or in plasmas with relativistic beams (P j m j c, where P j is the drift momentum). Using the relativistic expression for Ω j in Equation (2.4) shows that we can write for the pole of the function under the integral sign in Equation (3.1) is the Lorentz factor. We define the dimensionless parallel momentump ≡ p /m j c. The dimensionless parallel momentum associated with the relativistic pole is given bȳ We apply the technique proposed by Lerche (1967) to transform Equation (2.9) from the (p ⊥ , p ) coordinate system to the (Γ,p ) coordinate system (see also Swanson 2002;Lazar & Schlickeiser 2006;López et al. 2014López et al. , 2016. This transformation yields , and the Bessel functions are evaluated as J n ≡ J n (zp ⊥ ). Whenever ALPS performs a relativistic calculation and the code automatically transforms from (p ⊥ , p ) to (Γ,p ) coordinates and applies the polyharmonic spline algorithm described in Appendix C to create an equally spaced and homogeneous grid in (Γ,p ) coordinates. In this coordinate system, we perform the integration near poles and the analytic continuation in the same way as described in Sects. 3.1 and 3.2, but using the relativistic parallel momentum associated with the pole from Equation (

Test Cases and Results
In this section, we compare ALPS with known reference cases based on either our own or previously published results. For the calculations shown on the left, we keep k ⊥ dp = 10 −3 constant and scan through k . For the calculations shown on the right, we keep k dp = 10 −3 constant and scan through k ⊥ . The A/IC mode in quasi-perpendicular propagation corresponds to the kinetic Alfvén wave (KAW) at k ⊥ dp 1/ β 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 ωr of the frequency and its imaginary part γ.

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 2015) 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 0j 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 ω r and its imaginary part γ as functions of the parallel and perpendicular wavenumbers, respectively. We use β ⊥j = β j = 1 for both protons and electrons, and v A /c = 10 −4 , where β ⊥j ≡ w 2 ⊥j /v 2 A and β j ≡ w 2 j /v 2 A . We normalise all frequencies in units of the proton cyclotron frequency Ω 0p and all length scales in units of the proton skin depth d p ≡ v A /Ω 0p . The momentumspace resolution for the ALPS calculation in the quasi-parallel limit is n ⊥ = 320, n = 640, P max, p = 8m p v A , and P max, e = 0.19m p v A . In the quasi-perpendicular limit, we use n ⊥ = 240, n = 480, P max, p = 6m p v A , and P max, e = 0.14m p v A . In both cases, we set P max, j = P max,⊥j , J max = 10 −45 , n I = 5, n p = 100, and T lim = 0.01. We study the accuracy of the results depending on the resolution in Appendix 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.
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 |det D|, where D is the tensor from Comparison of dispersion maps from PLUME (left) and ALPS (right) for k ⊥ dp = k dp = 10 −3 . The lines show isocontours of constant lg |det D|. Minima in these maps correspond to solutions to the hot-plasma dispersion relation.
Equation (2.12), in the ω r -γ plane. They are a useful tool to find the initial guesses for ω r and γ for the Newton-secant root-finding search. Although the calculation of a dispersion map still requires the calculation of all χ j , it does not entail the application of the Newton-secant root-finding algorithm. Solutions to the hot-plasma dispersion relation appear as minima in these diagrams. We use a Maxwellian plasma model with β ⊥j = β j = 1 for both protons and electrons, k ⊥ d p = k d p = 10 −3 , and v A /c = 10 −4 . For the ALPS calculation, we use n ⊥ = 240, n = 480, P max, p = P max,⊥p = 6m p v A , P max, e = P max,⊥e = 0.14m p v A , J max = 10 −45 , n I = 5, n p = 100, and T lim = 0.01. Both the ALPS and the PLUME calculations reveal seven solutions to the dispersion relation. We note that the point ω r = γ = 0 is a maximum and does not represent a solution to the dispersion relation. The solutions at ω r = ±10 −3 Ω 0p and γ = −2.3 × 10 −10 Ω 0p are the forward and backward propagating A/IC waves. The solutions at ω r = ±2 × 10 −3 Ω 0p and γ = −5.4 × 10 −5 Ω 0p are the forward and backward propagating FM/W waves. The solutions at ω r = ±1.2 × 10 −3 Ω 0p and γ = −7.3 × 10 −4 Ω 0p are the forward and backward propagating slow waves (ion-acoustic waves). Lastly, the solution at ω r = 0 and γ = −7.2 × 10 −4 Ω 0p is the non-propagating slow mode, which is sometimes denoted 'entropy mode', (Verscharen et al. 2016(Verscharen et al. , 2017. 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.

Anisotropic Bi-Maxwellian
Distributions PLUME, like most other standard hot-plasma dispersion-relation solvers, also 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 1961;Davidson & Ogden 1975;Yoon et al. 2010), and as an example for a non-propagating instability, we calculate the dispersion relation for the mirrormode instability (Rudakov & Sagdeev 1961;Tajiri 1967;Southwood & Kivelson 1993). The thresholds for both of these instabilities fulfil T ⊥p > T p . For this demonstration, we  Figure 3. Comparison of dispersion relations for the A/IC instability (left) and the mirror-mode instability (right) from PLUME and ALPS. We use T ⊥p /T p = 3. For the calculation of the A/IC instability, we keep k ⊥ dp = 10 −3 constant and scan through k . For the calculation of the mirror-mode instability, we keep θ = 75 • constant and scan through |k|.
use PLUME to calculate ω r and γ as functions of the wavenumber in a plasma with bi-Maxwellian protons and Maxwellian electrons using β p = β e = β ⊥e = 1, T ⊥p /T p = 3, and v 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 ⊥ = 320, n = 640, and P max, p = 8m p v A , P max,⊥p = 13.9m p v A , P max, e = P max,⊥e = 0.19m p v A , J max = 10 −45 , n I = 5, n p = 100, and T lim = 0.01. We study the accuracy of these results depending on the resolution in Appendix 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.

Anisotropic κ-Distributions
Astfalk et al. (2015) developed the code DSHARK to calculate dispersion relations in plasmas with bi-κ-distributions. As one example, these authors discuss the FM/W instability in an anisotropic electron-proton plasma with κ p = κ e = 8, β p = 2, β e = 4, T ⊥p /T p = 0.4, and T ⊥e /T e = 0.5 (see Figure 1 from Astfalk et al. 2015). The angle between k and B 0 is constant for this calculation and set to θ = 0.001 • . We use DSHARK to reproduce this test case and set up κ-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 ⊥ = 400, n = 800, P max, p = 10m p v A , P max,⊥p = 6.32m p v A , P max, e = 0.33m p v A , P max,⊥e = 0.23m p v A , J max = 10 −45 , n I = 5, n p = 500, T lim = 0.01, and v A /c = 10

Relativistic Jüttner Distributions
As one example for a dispersion relation in a relativistic plasma, we reproduce the results by López et al. (2014) for an electron-positron pair plasma with a Jüttner distribution using v A /c = 1, m p = m e , β p = β e = (0.2, 0.4, 1.0) and T ⊥j = T 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 ⊥ d p = 10 −3 . We use n ⊥ = 30, n = 60, P max = 5m p v A , J max = 10 −45 , n I = 5, n p = 300, and T lim = 0.01. Our interpolation method transforms the (p ⊥ , p ) grid to the (Γ,p ) grid with n Γ = 500 and np = 500 steps in Γ andp , respectively. We show the results in Figure 5. López et al. (2014) 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. (2014) and ALPS is only visible in the real part of the frequency at the large-k /low-ω r end of the A/IC branches.

Discussion and Conclusions
ALPS solves the relativistic and non-relativistic hot-plasma dispersion relations in a plasma with arbitrary background distribution functions. We have benchmarked ALPS against existing codes by comparing dispersion relations for waves and instabilities in Maxwellian, bi-Maxwellian, κ-distributed, and relativistic Jüttner-distributed plasmas. In all cases, we find that ALPS agrees well with existing codes. This finding encourages us to apply ALPS to yet unexplored plasma environments in future work. An important application of ALPS will be the analysis of distribution functions measured by spacecraft in the solar wind. ALPS includes the necessary numerical framework to preprocess and format the spacecraft data so that they can serve as f 0 tables for direct input (see Appendix C). Especially, the upcoming missions Solar Orbiter and Parker Solar Probe will deliver plasma measurements with unprecedented energy and time resolution in the solar wind that will serve as the ideal input for ALPS. The vast majority of previous kinetic studies of waves and instabilities relied on bi-Maxwellian fits to the observed distribution functions and the use of a standard bi-Maxwellian code to solve the hot-plasma dispersion relation such as WHAMP, NHDS, or PLUME. Our approach allows us, however, to relax the bi-Maxwellian assumption and to analyse the plasma behaviour more realistically. Future comparisons of the results from standard codes such as PLUME with the results from ALPS will help to evaluate the quality of the previous bi-Maxwellian approaches and to refine our understanding of the role of instabilities in collisionless plasmas based on the actual distribution functions. For instance, our knowledge of the realistic value of certain instability thresholds is still very limited. Some in-situ observations of kinetic plasma features in the solar wind lie above the thresholds of kinetic instabilities when calculated based on bi-Maxwellian background distributions (see, for example, Isenberg 2012). The general conjecture is, however, that the plasma is limited by the lowest instability threshold. A more realistic calculation based on the actual distribution functions may resolve this discrepancy. This concept applies, for example, to anisotropy-driven instabilities such as the A/IC instability (Hellinger et al. 2006;Bale et al. 2009;Maruca et al. 2012) or beam-driven instabilities such as the FM/W instability (Reisenfeld et al. 2001;Verscharen et al. 2013a). Also nonthermal electron configurations, which are known to carry a significant heat flux into the solar wind, require a non-bi-Maxwellian representation for the determination of the relevant instabilities that limit their heat flux (Feldman et al. 1975;Pilipp et al. 1987a;Pulupa et al. 2011;Salem et al. 2013). Another field of application of ALPS is the study of highly non-thermal plasma configurations related to reconnection events (Phan et al. 2006;Gosling 2007;Gosling et al. 2007;Egedal et al. 2012Egedal et al. , 2013. We also emphasise the applicability of ALPS for the determination of dispersion relations using distributions from numerical plasma simulations. Particle-in-cell or Eulerian plasma codes generate data directly suitable as f 0 tables for ALPS. Some of these numerical simulations use (realistically or artificially) relativistic plasma conditions. Therefore, ALPS's ability to include relativistic effects will be very useful for the study of the wave properties and the stability of simulated plasmas.
Our resolution studies in Appendix A offer some insight into the necessary resolution of the f 0 tables for a reliable determination of the plasma dispersion relation. In the shown applications, a minimum resolution of about n ⊥ = 40 and n = 80 has proven to be necessary for a good agreement between ALPS and the test results for (bi-)Maxwellian distributions. In a future extension of ALPS, we will include Nyquist's method to automatically determine the stability of directly observed distribution functions (Klein et al. 2017).
We appreciate helpful comments and contributions from Sergei Markovskii and Thomas Brackett. The ALPS collaboration appreciates support from NASA grant NNX16AG81G. We present more details about the numerics on the website www.alps.space. The ALPS source code will be made publicly available on this website after our initial science phase. Computations were performed on Trillian, a Cray XE6m-200 supercomputer at UNH supported by the NSF MRI program under grant PHY-1229408. D.V. was supported by the STFC Ernest Rutherford Fellowship ST/P003826/1. B.D.G.C. was supported in part by NASA grants NNX15AI80G and NNX17AI18G and NSF grant PHY-1500041.

Appendix A. Resolution Studies
In order to understand the required resolution of the f 0 tables for calculations with ALPS, we compare results from PLUME with results from ALPS for the same plasma parameters using different resolutions in this appendix. For all calculations, we use β p = β e = 1, T ⊥e /T e = 1, J max = 10 −45 , n I = 5, n p = 100, and T lim = 0.01. We use P max, p as a free parameter and set n = 2n ⊥ , P max, e = P max,⊥e = P max, p m e /m p , and P max,⊥p = P max, p T ⊥p /T p . We define the resolution in momentum as ∆w j ≡ P max, j /(n m j w j ) and the resolution in frequency as ∆ω r /ω r ≡ |ω r,ALPS − ω r,PLUME | /ω r,PLUME , where ω r,ALPS is the solution from ALPS, and ω r,PLUME is the solution from PLUME. For κ-distributions, the appropriate resolution depends on both β j and κ. Instead of giving general guidelines for the resolution, we, therefore, recommend case-by-case convergence studies when calculating dispersion relations in plasmas with κ-distributions.

A.1. Maxwellian Distributions
In Figure 6, we show a resolution study for the A/IC wave in quasi-parallel propagation in an isotropic Maxwellian plasma. This figure complements our solutions shown in Figure 1. The four panels represent different values of P max, j . In each panel, the diagram at the top compares the real part of the frequency from five ALPS calculations with different ∆w j to the Maxwellian solutions from PLUME. The diagram at the bottom compares the ratio between ω r from the five ALPS calculations and ω r from PLUME. In this parameter range, P max, p = 8m p w p with a resolution finer than ∆w j = 0.1 leads to a very good agreement with the PLUME solutions for ω r . For wavenumbers below 1/d p , a lower value of P max, p is sufficient. Figure 7 shows the same as Figure 6, but giving the . Resolution study for the real part of the frequency for the A/IC-wave solution in quasi-parallel propagation. We keep k ⊥ dp = 10 −3 constant and scan through k .
imaginary part of the frequency instead of its real part. This figure confirms our finding regarding the optimal resolution. Figures 8 and 9 show the same as Figures 6 and 7, but for quasi-perpendicular propagation instead of quasi-parallel propagation. The required resolution is lower in the quasi-perpendicular case than in the quasi-parallel case. The solutions with P max, p = 4m p w p and ∆w j 0.1 lead to a very good agreement between the ALPS and PLUME solutions.

A.2. Anisotropic Bi-Maxwellian Distributions
In addition to our Maxwellian test, we study the dependence of the ALPS solutions on the resolutions for the bi-Maxwellian case with T ⊥p /T p = 3 as shown in Figure 3. Figure 10 compares ALPS solutions for the A/IC instability in quasi-parallel propagation for different values of P max, p and ∆w j with the solutions from PLUME for the real part of the frequency. Figure 11 compares ALPS and PLUME solutions for the imaginary part of the frequency. The solutions with P max, p = 4m p w p and ∆w j 0.05 lead to a good agreement between ALPS and PLUME in both ω r and γ.
In Figure 12, we study the dependence of the solutions on the resolution for the mirrormode instability with the same parameters as in Figure 3. The correct solution of the mirror-mode instability has ω r = 0; however, the ALPS solutions have finite values ω r = 0. The value of ω r decreases with increasing ∆w j . As Southwood & Kivelson (1993) point out, the mirror-mode instability is strongly influenced by particles with p ≈ 0. The error in frequency ∆ω r is determined by the resolution of the momentum grid around Resolution study for the imaginary part of the frequency for the A/IC-wave solution in quasi-parallel propagation. We keep k ⊥ dp = 10 −3 constant and scan through k . p = 0, where ∆ω r ∼ k w j ∆w j . Figure 13 shows the comparison of the imaginary part of the mirror-mode solutions. Like in the case of the A/IC instability, a resolution with P max, p = 4m p w p and ∆w j 0.05 leads to a good agreement between ALPS and PLUME.

Appendix B. Levenberg-Marquardt Fit
For the hybrid analytic continuation, ALPS fits the f 0 table with a combination of pre-described algebraic expressions as described in Section 3.2. We employ a Levenberg-Marquardt algorithm (Levenberg 1944;Marquardt 1963) to fit the distribution functions in p with a superposition of an arbitrary number of Maxwellian distributions, κdistributions, and Jüttner distributions. The user can freely choose the number of fits and their superposition. We evaluate different fit parameters for each given value of p ⊥ . We define the Maxwellian fitting function as where u k are the fit parameters, y is a constant user-defined parameter, andp ⊥ andp are the normalised perpendicular and parallel momenta. The parameter y compensates the otherwise strong p ⊥ -dependence of u 1 , making the fit more reliable. It is constant for all p ⊥ . We choose this expression rather than a fit in p ⊥ since it provides a greater flexibility in the p ⊥ -domain compared to a two-dimensional fit in p ⊥ and p . The best choice for y is β ⊥j m p /m j . The standard normalisation in ALPS  Figure 8. Resolution study for the real part of the frequency for the A/IC-wave solution in quasi-perpendicular propagation. We keep k dp = 10 −3 constant and scan through k ⊥ . andp = p /m p v A . In cases with κ-distributed plasma components, we use .

(B 2)
In cases with Jüttner-distributed plasma components, we use In this case, the best choice for y is 2c 2 /w 2 j . These fitting relations are easily extendable by the user to cover more general functions as needed.
We denote the discretised f 0 table of species j at constant p ⊥ asf i,j (p ,i ), the discrete steps inp asp ,i , the vector of all fit parameters as u, and the sum of all fit functions as F (p ). In the Jüttner-distributed cases, the coordinates are replaced with Γ andp accordingly. We define the residuals as s i ≡f i,j − F (p ,i ) and define C ≡ i s 2 i . We denote the Jacobian off j with respect to u as J. We use a superposition of analytical expressions for the Jacobian based on the given form of F (p ).
The Levenberg-Marquardt algorithm uses an iterative step to update u of the form where λ is a user-defined scalar. For the matrix inversion in Equation (B 4), we use the LU-factorisation. Then we calculate the residuals s new based on u new and determine C new = |s new | 2 . If C new C, we set u to u new , reduce λ by a constant factor λ f (userdefined, standard value is 10), and repeat the procedure. If C new > C, we discard u new , increase λ by the constant factor λ f , and repeat the procedure. In this way, we iteratively . Resolution study for the imaginary part of the frequency for the A/IC-wave solution in quasi-perpendicular propagation. We keep k dp = 10 −3 constant and scan through k ⊥ .
determine the fit parameters u until the fit converges (i.e., C with a user-defined ), or until the number of iterations reaches a user-defined maximum value. ALPS writes the fitted distribution into a separate output file so that a direct comparison with the original input distribution is possible.

Appendix C. The Smoothed Thin-Plate Spline Interpolation
Spacecraft or other plasma data are typically not available on a dense Cartesian grid like the grid required for an f 0 table in ALPS. Therefore, our code includes an interpolation algorithm that fills gaps between data points. ALPS uses the same interpolation algorithm to create an equidistant grid in (Γ,p ) space after the coordinate transformation in cases with relativistic poles. We use a polyharmonic spline interpolation with the radial basis function of a thin-plate spline with smoothing (Powell 1994;Donato & Belongie 2002). For each species, we begin with the "coarse" distribution function f c,µ which is given by n c data points (index µ = 1 . . . n c ) with the associated coarse momentum coordinatesp ⊥c,µ andp c,µ . The set (f c,µ ,p ⊥c,µ ,p c,µ ) forms one data point. The coarse grid is typically not equally distributed in momentum space.
For each species, the "fine" grid of momentum coordinates is given byp ⊥,i,k andp ,i,k with i = 1 . . . n ⊥ and k = 1 . . . n (and correspondingly in the coordinates Γ andp for cases with relativistic poles). The fine grid corresponds to the actual f 0 table to be used as input in ALPS. The goal of our interpolation is to find the value of the distribution functionf i,k on all grid points (i, k). We define the vectors w = (w 1 , . . . , w nc ), c = . Resolution study for the real part of the frequency for the A/IC-instability solution in quasi-parallel propagation. We use a bi-Maxwellian plasma with T ⊥p /T p = 3. We keep k ⊥ dp = 10 −3 constant and scan through k .
(c 1 , c 2 , c 3 ),f c = (f c,1 , . . . ,f c,nc ), and 0 = (0, 0, 0). We furthermore define the matrix K µ,ν = r 2 log(r) if r 1 r log(r r ) if r < 1, (C 1) where r ≡ (p ⊥c,µ −p ⊥c,ν ) 2 + (p c,µ −p c,ν ) 2 . We also define the (n c × 3) matrix P. Its µth row is given by (1,p ⊥c,µ ,p c,µ ). The thin-plate spline interpolation requires to solve the nonhomogeneous linear system of equations for the vectors w and c. α is a user-defined smoothing parameter (α = 0 forces the fine grid to run through all points of the coarse grid), and 1 is the (n c × n c ) unit matrix. The interpolation is then given bŷ f i,k = c 1 + c 2p⊥,i,k + c 3p ,i,k + nc µ=1 w R µ i,k , (C 3) where R µ i,k ≡ (p ⊥,i,k −p ⊥c,µ ) 2 + (p ,i,k −p c,µ ) 2 . The numerically expensive part of the interpolation is the solution of Equation (C 2). Since K 11 = 0, a direct LU fac- . Resolution study for the imaginary part of the frequency for the A/IC-instability solution in quasi-parallel propagation. We use a bi-Maxwellian plasma with T ⊥p /T p = 3. We keep k ⊥ dp = 10 −3 constant and scan through k .  Figure 13. Resolution study for the imaginary part of the frequency for the mirror-mode-instability solution. We use a bi-Maxwellian plasma with T ⊥p /T p = 3. We keep θ = 75 • constant and scan through |k|.
torisation is not possible. Therefore, we apply a LU-factorisation algorithm with partial pivoting through row permutations until K 11 = 0.