A linear dispersion relation for the hybrid kinetic-ion/fluid-electron model of plasma physics

A dispersion relation for a commonly used hybrid model of plasma physics is developed, which combines fully kinetic ions and a massless-electron fluid description. Although this model and variations of it have been used to describe plasma phenomena for about 40 years, to date there exists no general dispersion relation to describe the linear wave physics contained in the model. Previous efforts along these lines are extended here to retain arbitrary wave propagation angles, temperature anisotropy effects, as well as additional terms in the generalized Ohm's law which determines the electric field. A numerical solver for the dispersion relation is developed, and linear wave physics is benchmarked against solutions of a full Vlasov-Maxwell dispersion relation solver. This work opens the door to a more accurate interpretation of existing and future wave and turbulence simulations using this type of hybrid model.


Introduction
Phenomena such as plasma instabilities, magnetic reconnection, turbulent transport and dissipation of energy have been studied for several decades through experiment and observation, and in a variety of systems ranging from laboratory plasmas like linear devices and tokamaks to natural plasmas like the solar corona, solar wind, or planetary magnetospheres (e.g., Refs. [1][2][3][4][5]). The plasmas permeating these systems are in many cases sufficiently collisionless to support a host of kinetic effects, with strong influence on the evolution of those systems. A comprehensive theoretical description of such a plasma often proves to be very complex, as analytical solutions of the kinetic equations are usually constrained to simplified cases or linearized problems. Numerical methods must then be employed in order to model more realistic situations.
However, because of the intrinsic separation between ion and electron spatiotemporal scales in a plasma, in many cases numerically simulating the full Vlasov-Maxwell system of equations for all involved species is also not feasible due to the significant computational expense of evolving a six-dimensional system (3 spatial, 3 velocity-space dimensions) across the diverse timescales that are involved. In practice, many simplifications such as reduced spatial dimensionality, an artificially reduced ion/electron mass ratio, or analytically reduced models such as gyrokinetics [6][7][8] are employed to make such investigations tractable.
An alternative and very common approach along these lines is to introduce a hybrid model that applies the fully kinetic treatment only to the ions, while describing the electrons within the framework of a simpler fluid model. Historically, such models have been used for about 40 years [9][10][11][12][13][14][15], but surprisingly, to our knowledge there exists no general dispersion relation to describe the linear wave physics within such a reduced model. Limiting cases for perpendicular propagation [16] and for parametric instabilities of Alfvén waves propagating parallel to the background magnetic field [17] have been derived, but the general case of oblique propagation has not been treated in the literature. Existing simulation codes based on this model (see, e.g., Refs. [18][19][20][21][22][23] for recent efforts) are often benchmarked against frequencies (more rarely against growth or damping rates) obtained from analytical or full Vlasov-Maxwell dispersion relations.
In the present paper, we intend to fill this gap by deriving a dispersion relation for the hybrid kinetic ion, massless electron model, variations of which are being actively used in a number of simulation codes throughout the space and astrophysics community. The existence of such a dispersion relation will enable a more thorough and clear-cut validation of the hybrid model itself, and also of simulations making use of that model. For this purpose, we have developed the numerical dispersion solver "HYDROS", which is publicly available [24] and has already been used to help interpret nonlinear simulation results in a recent publication [25].
This paper is structured as follows: In Sec. 2, we describe briefly the set of equations that will form the basis for our dispersion relation. In the main part of the paper, Sec. 3 is devoted to the derivation of the dispersion relation, while in Sec. 4 we describe the implementation of the "HYDROS" solver, together with a variety of benchmark tests against the fully kinetic code DSHARK [26]. Sec. 5, finally, provides a summary of this work.

The hybrid-kinetic system of equations
In this section, we introduce the hybrid kinetic-ion/fluid-electron ("hybrid-kinetic" in the following) equations that will be used throughout this paper. For our present purposes, we stick to a nonrelativistic, low-frequency (i.e. ω Ω ce ) version of the equations which assumes massless electrons and retains only a singly charged ion species such that n i = n e . The system of equations then consists of the ion Vlasov equation, and an Ohm's law which determines the electric field, with n i u i = vf i d 3 v, the resistivity η, and the electron pressure gradient ∇P e = C∇n γ e . The electromagnetic fields are further constrained by Faraday's law and the nonrelativistic version of Ampere's law These equations contain the full ion kinetic physics including wave-particle interactions such as Landau and transit-time damping, as well as cyclotron resonances. The ion background distribution is assumed to be bi-Maxwellian, enabling the occurrence of the basic ion-anisotropy driven instabilities such as the firehose and mirror modes. Electrons, on the other hand, appear only as a neutralizing, massless background species (where C = n 1−γ e T 0e , and the choice γ = 1 or 5/3 results in an isothermal or adiabatic electron model for the electron pressure gradient term) and implicitly as the carriers of the current j. Interactions between waves and electrons are not retained, so that, e.g., electron Landau damping is absent.

Derivation of the hybrid-kinetic Vlasov-Maxwell dispersion relation
This section details the derivation of a dispersion relation which completely describes the linear wave physics contained in the system defined by Eqs. (1)-(4). Here, we will loosely follow the procedure employed before for perpendicularly propagating waves in a hybridkinetic model in Ref. [16], although from the outset we retain a general wave vector, as well as additional terms in Ohm's law and an anisotropic background distribution. We note that alternatively, it is possible to derive a hybrid-kinetic dispersion relation using the dielectric tensor method described, e.g., in Ref. [27]. Since the dielectric is additive in species, it is possible to combine the fully kinetic ion susceptibility with that of a fluid electron model, derived, e.g. in Ref. [28].
As a first step, we linearize Eqs. (1)-(4) with respect to a static, homogeneous, background, according to the rules f i = F 0i + δf i , n i = n 0i + δn i , B = B 0 z + δB, and E = δE. The ion density perturbation is defined as δn i = δf i d 3 v. Furthermore, we introduce a plane wave expansion for all perturbed quantities For the purposes of this paper, we will assume the background ion distribution to have the shape of a gyrotropic Maxwellian, retaining the effects of an anisotropic temperature. The background distribution can thus be written as Under such an assumption, we may choose the alignment of our coordinate system such that, without loss of generality, The divergence-free property of the magnetic field then reads k · δB = 0, so that the magnetic fluctuations in x and parallel (z) direction are related via and we need to describe only δB y and δB separately. Applying the linearization rules to Eqs. (1)-(3), setting n e = n i , and using Eq. (4) to replace the current yields Because we chose a static background, we may now write u i = vδf i d 3 v/n 0i . In Eq. (7), the term proportional to v × B 0 · ∇ v F 0i vanishes because of the gyrotropic background distribution. The general procedure now is to insert the plane wave expansion into the above equations, and then derive a system of equations for δn i , δB y and δB , using Eqs. (8) and (9) to eliminate any dependencies on the other perturbed quantities. In the following, we will use only the Fourier coefficients, omitting the k index and the tildes denoting the perturbed quantities.
As a first step, we may apply the Fourier expansion to Eq. (8) and insert the result into Faraday's law, Eq. (9). This allows us to eliminate the electric field δE and, after some algebra, to solve for Next, we expand the Vlasov equation, Eq. (7), in plane waves, yielding Note that the term proportional to v×δB·∇ v F 0i in Eq. (11) vanishes only in an isotropic Maxwellian plasma, and thus must be kept for the present purposes. Introducing cylindric coordinates v ⊥ cos θ, v ⊥ sin θ, v in velocity space, we can obtain the relations Inserting these, the linearized Ohm's law from Eq. (8), and Eq. (10) into Eq. (11), we may write Now, we define a function g (θ) = RHS [Eq. (12)] /Ω c , so that this equation reduces to the form where we set We may now obtain a solution for Eq. (13) by means of the integrating factor method, as described similarly also in Ref. [16]. The integrating factor is defined as with ν = ω − k v /Ω ci and κ = k ⊥ v ⊥ /Ω ci . The solution for the distribution function is then given by Here, we introduced the abbreviation to make the following equations more concise. Note that we chose the lower integration boundary as −σ∞, where we set σ = 1 when (ω) > 0 and σ = −1 when (ω) < 0. This method ensures that the integration does not contain diverging parts, but at the same time precludes convergence of the integral for marginally stable modes. In the numerical solution of the dispersion relation, we will thus take steps to ensure that the solver cannot move arbitrarily close to zero (see Sec. 4).
Using the above expression for the ion distribution function, we may write the perturbed density as the velocity space integral of Eq. (14), which constitutes the first equation of the matrix system we will use to determine the wave solutions of the hybrid-kinetic system. Two more equations may be obtained from the x and y component of Eq. (10). The former yields and from the latter we can calculate In order to obtain a solution for this system of equations, we need to solve a set of integrals where a ∈ {0, 1}. We choose b ∈ {0, 1, 2} to denote the three choices for f (θ) ∈ {1, v ⊥ cos θ, v ⊥ sin θ}, and c ∈ {0, 1, 2} to denote the three choices for h (θ ) ∈ {1, v ⊥ cos θ , v ⊥ sin θ }. Finally, the velocity space volume element is given by d 3 v = v ⊥ dv dv ⊥ dθ, and V is used to denote the velocity space itself. Our present system of equations requires 15 of the 18 combinations defined thus. The solutions of these integrals can be obtained in the same fashion as described in Ref. [16], and is detailed in Appendix 5. Once the integrals are solved, we can express Eqs.

Introducing the abbreviations
and simplifying also Eqs. (15)-(17), we finally obtain the matrix elements The dispersion relation of the hybrid-kinetic model is then obtained by setting and solving for the complex frequencies that fulfill this equation.

Numerical implementation
The dispersion relation derived in Sec. 3 has been implemented for numerical solution using the Python/NumPy/SciPy framework, whose mathematical library provides all of the necessary functionality. Both the plasma dispersion function and the modified (and exponentially scaled) Bessel functions which appear when solving the integrals of Eq. (18) are provided by SciPy [29] via the specfun library [30], and can thus be readily used. In order to find the zeros of the dispersion relation, we employ the root finding methods provided by the SciPy library, with a Levenberg-Marquardt algorithm [31,32] set as the default method.
To improve the convergence speed of the root finding algorithm, several predictive algorithms are available to anticipate the evolution of a root during a parameter scan, or to direct the solver closer towards a particular (e.g. less damped) solution when several roots exist close together. These algorithms are a) using the old position of the root as a starting point for the new iteration, b) quadratic extrapolation using the last two positions of the root, c) modified quadratic extrapolation. The last method consists of a quadratic extrapolation with a subsequent modification applied to the predicted frequency or damping rate and can, e.g., be used for cases with degenerate modes. For all cases shown here, the dispersion relations were sufficiently clear that method b) could be employed without additional modifications. However, we found it very useful to incorporate diagnostics to ensure that the solver keeps following a particular solution, e.g., by comparing complex frequencies, field amplitudes, and cross phases between the fields.
As mentioned in Sec. 3, the integrals solved here do not converge if the imaginary part of the frequency is zero. For that reason, we introduce a lower boundary γ min = 10 −13 for the absolute value of the imaginary part. If the solver converges to a number smaller than γ min , the starting point for the next solution will be reset to γ min . This prevents a runaway of the solver towards ever smaller imaginary parts that has been observed otherwise.

Benchmark of HYDROS against other kinetic solvers
In this section, we aim to verify both the derivation of Sec. 3 and the numerical implementation of the HYDROS code by comparing its results to solutions of the fully kinetic dispersion solver DSHARK [26]. For this exercise, we choose several example waves and instabilities that are commonly encountered in many systems of interest (e.g., the solar wind and planetary magnetospheres), such as the fast and slow magnetosonic mode, the kinetic Alfvén wave (KAW), the oblique and parallel firehose mode and the mirror mode. For all studies in this work, we set the polytropic coefficient γ = 1 to obtain an isothermal electron description, and we choose the resistivity η = 0. This section will be subdivided by the propagation angle of the examined modes. We note that neither DSHARK nor HYDROS are able to treat exact parallel and perpendicular propagation, as these cases lead to divisions by zero in the general expressions for the dispersion relation. "Parallel" and "perpendicular" thus refer to an angle very close to 0 degrees or 90 degrees, respectively. In DSHARK, we furthermore set v A /c = 0.01 throughout (this ratio is undefined in HYDROS, as the displacement current is neglected).
In this section, wavenumbers will be plotted in units of kd i , where d i = c/ω pi = m i c 2 /4πn i e 2 , and frequencies will be given in terms of the ion cyclotron frequency Ω ci = eB 0 /m i c. The plasma beta parameter for a species s will be defined as β s = 8πn s T s /B 2 0 .

Parallel propagation
Parallel firehose instability. As the first example, we compare growth rates and frequencies of the parallel firehose instability, which occurs when the parallel thermal pressure significantly exceeds the perpendicular thermal pressure. Here, we choose the parameters β p = 4, β p⊥ = β e = 1, and we set the propagation angle to 10 −4 degrees. The results for these parameters from both the fully kinetic and the hybrid-kinetic solver are shown in Figs. 1a and 1b. Both in growth rates and frequencies we find very good agreement between the two solvers.
Ion cyclotron instability. Inverting the anisotropy used for the parallel firehose instability to β p = β e = 1, β p⊥ = 4 while keeping the propagation angle parallel to the background field yields a plasma that is susceptible to the ion cyclotron instability. Scanning over the wavenumber (see Figs. 2a and 2b) then reveals a growing wave with a frequency close to the ion cyclotron frequency, and a growth rate peak close to kd i = 0.7. As before, DSHARK and HYDROS agree very well on both growth rates and frequencies.
L-mode. Next, we choose an isotropic ion distribution and examine the results for the L-mode solution (Alfvén/ion-cyclotron wave). Here, we select two different values of the plasma β, namely β i = β e = 10 −2 and β i = β e = 1. The relevance of the β value in this case is its influence on the effect of the ion cyclotron resonance -for low β, one indeed finds that the wave frequency asymptotes close to Ω ci , whereas it converges towards a lower frequency in the finite-β case, accompanied by a significantly stronger damping. Once again, both DSHARK and HYDROS agree very well on both this effect and on the dispersion curves for the two cases, which are plotted in Figs. 3a and 3b. Furthermore, for low β, very good agreement with an analytical fluid dispersion relation [19] is found.
R-mode. The R-mode (or whistler) is characterized by its dispersive behavior, with frequencies (for parallel propagation) proportional to k v A for k d i ≤ 1, but proportional to k v A 2 at higher wavenumber. In Figs. 4a and 4b the damping rates and frequencies for this wave are presented, for β i = β i⊥ = β e = 1, along with an analytical dispersion relation [19]. A small frequency deviation between DSHARK and HYDROS towards higher wavenumbers can be observed, which follows from neglecting electron    Figure 4: R-mode damping rates (a) and frequencies (b) for HYDROS and DSHARK, using β i = β i⊥ = β e = 1. For comparison, an analytical R-mode frequency formula [19] is plotted for the limit of zero electron mass and low β (blue).
inertia effects in the hybrid-kinetic model. The latter thus shows a perfect k 2 scaling towards higher wavenumbers, whereas the fully kinetic model eventually rolls over and asymptotes towards close to Ω ce (not shown). This effect is also kept in the analytical fluid dispersion relation [19] plotted as well in Fig. 4b. In the damping rates of Fig. 4a, a slight discrepancy can be discerned around kd i ≈ 1, which may be attributed to the effects of electron Landau or Barnes damping [33], which is not retained within the hybrid-kinetic model. The electron mass may be artificially reduced towards the limit taken in our hybrid derivation, improving the agreement between the two curves (not shown). Ion acoustic waves. Finally, we examine the dispersion relation of the ion acoustic mode. Although the wavevector of this mode can be arbitrarily oriented with respect to the background field, we choose to examine it for parallel propagation here to enable a cleaner comparison to analytical theory, avoiding visible cyclotron damping effects. Figs. 5a and 5b present the code results, as well as analytical formulae for damping rates and frequencies [19]. Note that the latter apply for T i /T e 1, but a relatively large value of 0.04 has been chosen for this ratio to achieve well-measurable damping rates. Thus, although the results from both numerical solvers agree well with other, they deviate somewhat from the analytical damping rates. In addition, for this wave a very small electron mass m e = 10 −12 m i had to be chosen in DSHARK to avoid electron Landau damping, which would otherwise dominate here.

Perpendicular propagation
Ion Bernstein modes. Next, we focus on purely perpendicularly propagating ion Bernstein modes, for β i = β i⊥ = β e = 1. For the first ion Bernstein mode close to ω = Ω ci , we find that its peak frequency deviates from the DSHARK result by about 5% (see Fig. 6a). This discrepancy, however, is resolved when approximating, in the kinetic dispersion solver, the massless electron limit that enters the HYDROS equations. Good numerical agreement is obtained for m e /m p 10 −12 , and for Figs. 6a and 6b a value of m e = 10 −15 m p has been used. With this setup, excellent agreement is obtained for the five first ion Bernstein modes, as demonstrated in Fig. 6b. The damping rates obtained for these modes are numerical zeros (Landau and transit time damping are negligible because k ≈ 0) and are thus not plotted.

Oblique propagation
In this subsection, we explore the general case of oblique propagation, in which various processes, e.g. resonant wave-particle interactions that depend on perpendicular and parallel velocity, occur in combination, resulting in more complex mode behavior.
Oblique firehose instability. The first instability we examine is the oblique firehose instability, which results in a non-propagating (zero frequency) mode whose growth rate depends on the wavevector. Here, we choose β i = 10, β i⊥ = 26/3, β e = 1 and a propagation angle of 45 • . In Fig. 7a, we compare the latter between HYDROS and DSHARK, showing excellent agreement of the growth rates. Frequencies are not shown due to the non-propagating nature of this mode.
Mirror mode. Next, we examine the mirror mode, which responds to an anisotropy of the ion distribution that favors the perpendicular pressure. We set β i = 6, β i⊥ = 10, β e = 1, keeping the propagation angle of 45 • . The resulting growth rates are shown in Fig. 7b and, once more, show good agreement between the two codes. A close look, however, shows some minor discrepancy close to the growth rate peak. This is again a result of the electron physics missing from the hybrid model, and can be resolved when taking the electron towards zero in the fully kinetic code.
Kinetic Alfvén wave. The kinetic Alfvén wave is the kinetic equivalent of the MHD shear-Alfvén wave, and will be examined here for a propagation angle of 85 • , and β e = β p = β p⊥ = 1. We scan the wavenumbers from kd i = 0.1 up to kd i = 20, using first the real proton/electron mass ratio in DSHARK (squares in Figs. 8a and 8b). As is demonstrated in Fig. 8a, we find very good agreement in the frequencies up to about kd i ∼ 3. Beyond this wavenumber, however, the hybrid model KAW converges to a frequency of about 1.5Ω ci , whereas in the kinetic model it stays slightly below Ω ci .
Although barely visible in Fig. 8b, we find that the hybrid model underpredicts the fully kinetic KAW damping at ion scales (below kd i ∼ 1) by about 25%, which can be attributed to the missing electron Landau damping. Above kd i ∼ 1, this gap becomes even more significant, but closes again as the ion cyclotron damping becomes dominant at about kd i ∼ 4. Reducing the electron mass in DSHARK to m e = 10 −12 m p , the hybrid results both for the frequency and the damping rates are recovered with very good accuracy (pentagons in Figs. 8a and 8b). Fast magnetosonic mode. Our next test concerns the fast magnetosonic mode, in the wavenumber range where its frequency approaches that of the ion cyclotron motion. For the chosen propagation angle of 85 • , this occurs roughly at a wavenumber of kd i ≈ 0.6. As before, we use β i = β i⊥ = β e = 1. For these parameters, the fast mode has a left-handed polarization, such that it resonates with the ion gyration and is confined to frequencies below the ion cyclotron frequency. We present the comparison of the resulting frequencies in Fig. 9a, and the corresponding damping rates in Fig. 9b.
In the frequency plot, good agreement is observed, although at low wavenumbers there is a slight shift between the curves of the hybrid and the fully kinetic model until the ion cyclotron frequency is approached. In the same wavenumber range, the damping rate plot exhibits more severe disagreement: for proton/electron mass ratio (squares) the fully kinetic model yields significant finite damping rates, whereas the hybrid model predicts undamped waves. In this case, the observed wave damping is caused by electron transit time damping [33], and proves to be rather resilient when changing the electron mass -in order to obtain good agreement with the hybrid model, as shown in Fig. 9b, the electron mass had to be reduced to m e = m p /10 6 . This modification, at the same time, removes the frequency shift observed before.
Slow magnetosonic mode. Finally, we study the dispersion relation of the slow magnetosonic mode as we transition from the MHD spatial scales into the kinetic range. As before, we use β i = β i⊥ = β e = 1, and a propagation angle of 85 • . We start the scans from kd i = 0.01, and both solvers are able to track this heavily damped mode down to kd i 10, where its damping rate exceeds the frequency by about a factor 10. In this case, no adjustment of the electron mass is required, and the HYDROS results agree very well with those of DSHARK, as is demonstrated in Figs. 10a and 10b.

Conclusions
In this work, we have derived a dispersion relation for a widely used hybrid kineticion/fluid-electron model of plasma physics, which comprises the full physics content of a kinetic ion description, while retaining only a simple massless electron fluid model. This dispersion relation is valid for arbitrary propagation angle, and retains basic anisotropy effects by allowing for a bi-Maxwellian background ion distribution, while focusing on the description of a single ion species.
We have described the implementation of the dispersion solver HYDROS, which utilizes the facilities of the Python/SciPy language to solve the hybrid dispersion relation numerically. The correctness of our derivation and implementation has been verified against the fully kinetic code DSHARK, for a variety of waves and instabilities that propagate parallel, perpendicularly, or obliquely with respect to the background magnetic field.
In all examined cases, the hybrid-kinetic dispersion solver showed the expected behavior and faithfully reproduced the ion physics contained also in a fully kinetic model. In some instances, it was necessary to numerically approach the massless electron limit in the fully kinetic code in order to obtain good agreement, indicating the impact of electron physics that is not captured in the present version of the hybrid-kinetic model. Although it was always possible in principle to recover the hybrid results by simply using very light electrons in DSHARK, the use of a dedicated hybrid solver proved to be numerically advantageous, since the occurrence of diverging Bessel arguments is avoided. Furthermore, HYDROS is suitable to serve as a testbed for additional linear physics introduced through more sophisticated electron fluid models, which are not obtained by taking the small mass ratio limit of full kinetics.
We emphasize that in the present paper we made no attempt to examine the quality or applicability of the hybrid model itself for real-world plasmas. Instead, our main purpose here was to demonstrate the correctness of both our derivation and its numerical implementation. In future work, we aim to fill this gap and give a more complete account of the physics fidelity of the hybrid model for both natural and laboratory plasmas, and its relative merit in comparison to other reduced models such as gyrokinetic theory. In particular, we expect such a study to have a bearing on turbulence studies for solar wind and magnetospheric plasmas, for which hybrid-kinetic simulation codes are a common workhorse. To this end, the HYDROS code has been made publicly available in a Github repository [24].

Appendix: Gyroangle integrals
In this appendix, we describe the solution method for the integrals that were introduced in Eq. (18) of Sec. 3. Their definition, written formally as results in 18 separate integrals. In order to solve the dispersion relation for the hybrid model equations as introduced in Sec. 2, we require only those 15 integrals (the ones defined by a = 0 ∧ c = 0 do not appear), which are explicitly written as All integrals can be calculated with the same approach, which will be demonstrated here for the first integral As a first step, we solve the innermost integral over θ by means of substitution. We set θ = θ − σφ and obtain Next, we write cos (θ − σφ) in exponential form and separate the integrand into a secular part and a periodic part h (φ). Then we replace the infinite φ integration boundaries by an infinite sum of 2π slices, i.e.
Applied to the integrand of I 001 this reads We introduce a shift φ = φ − 2πn: Next, we may perform the integration over θ, which yields Upon performing the φ integration, the geometric series Q (ν) as well as the factor σ cancel out and we obtain Now, it remains to perform the integrations in v and v ⊥ . We write and insert F 0i from Eq. (5) to obtain Here, we have introduced the plasma dispersion function [34] Z Finally, we perform the v ⊥ integration and, introducing a further abbreviation obtain the solution We have introduced the abbreviations Since the Γ m function and its derivative Γ m appear exclusively with the argument ζ, it is omitted in the following. The remaining integrals can be solved in the same fashion, and their final expressions read