Comparative study of gyrokinetic, hybrid-kinetic and fully kinetic wave physics for space plasmas

A set of numerical solvers for the linear dispersion relations of the gyrokinetic, the hybrid-kinetic, and the fully kinetic model is employed to study the physics of the kinetic Alfv\'en wave and the fast magnetosonic mode in these models. In particular, we focus on parameters that are relevant for solar wind oriented applications (using a homogeneous, isotropic background), which are characterized by wave propagation angles averaging close to 90{\deg}. It is found that the gyrokinetic model, while lacking high-frequency solutions and cyclotron effects, faithfully reproduces the fully kinetic Alfv\'en wave physics close to, and sometimes significantly beyond, the boundaries of its range of validity. The hybrid-kinetic model, on the other hand, is much more complete in terms of high-frequency waves, but owing to its simple electron model it is found to severely underpredict wave damping rates even on ion spatial scales across a large range of parameters, despite containing full kinetic ion physics.


Introduction
Plasmas both in nature or in the laboratory often exist in hot and/or dilute states, resulting in very long mean free paths of the charged particles. Under such conditions, traditional fluid approaches like magnetohydrodynamics are often still applicable to the large-scale evolution of the system, but they do not account for kinetic effects like wave-particle interactions. On the one hand, these may cause microinstabilities resulting in enhanced transport, and on the other hand, through effects like Landau and cyclotron damping they may determine how energy is dissipated at the end of the turbulent spectrum. The former issue, i.e. microinstabilities causing enhanced turbulent transport, is one of the key challenges in fusion research [1], whereas the latter question about the nature of energy dissipation has recently been termed "one of the outstanding open problems in space physics" [2].
In order to describe such phenomena, one would ideally solve the fully kinetic Vlasov-Maxwell system. Unfortunately, this is usually not feasible for a realistic system involving multiple scales from the global system size down to the Debye length scale. In order to circumvent this problem, often approximations are imposed, either by artificially reducing the dimensionality of such simulations, or by employing a reduced model that is computationally less intense. Such theories have been developed both for laboratory systems such as fusion plasmas, and for other systems like the solar wind, the corona, or astrophysical plasmas. While the gyrokinetic model [3] is a de-facto standard for turbulence modeling in core fusion plasmas-sometimes employing further optimizations such as a low mass ratio expansion for the electron dynamics [4,5]-, in space physics several different kinetic models are used by various groups to describe similar physics (see below for references). A recent publication [6] has sparked the "turbulent dissipation challenge" (similar in spirit to the GEM reconnection challenge [7]), which aims to compare the various reduced models to determine their capability of modeling various aspects of plasma physics that are relevant to the solar wind, but which should also prove instructive for neighboring fields. The present work aims to aid this effort by comparing, within a linear framework, the gyrokinetic model [8][9][10][11][12][13][14][15], a widely used hybrid-kinetic model [16][17][18][19][20][21][22][23][24][25][26][27][28], and a fully kinetic ion and electron description.
To date, such comparisons have been mainly performed between two models [8,29], but so far there exists, for instance, no direct comparison between the gyrokinetic and the hybrid-kinetic model for parameters relevant to the solar wind, a gap that we intend to fill with the present work. This paper is structured as follows: In Sec. 2, we give a brief account of the three models that will be compared in the present work. In the main part, Sec. 3, a detailed comparison between these three models is performed for two commonly analyzed waves, the kinetic Alfvén wave, and the fast magnetosonic mode/whistler. In Sec. 4, we summarize these results and provide some conclusions as well as an outlook.

The gyrokinetic system of equations
The gyrokinetic model has originated and been employed successfully for several decades in the fusion community, and its range of applications has only recently been extended to space plasmas. Gyrokinetics (GK) is an analytical limit of full kinetics which is intended for strongly magnetized plasmas, where the particle motion perpendicular to a strong magnetic guide field can be expanded in terms of a fast gyration around that field and a drift motion perpendicular to that field. Several ordering assumptions must be obeyed in order for this theory to be valid [3,8], i.e.
where is a small parameter. As a consequence of the above ordering relations, the wave propagation angles that can be described within GK are constrained to be almost perpendicular to the background field. In addition, wave frequencies formally need to be much smaller than the cyclotron frequency of the involved species. In the main part of this paper, we will analyze how much the wave physics is actually altered when one or more of these conditions is violated. The gyrokinetic dispersion relation for a homogeneous plasma with isotropic Maxwellian background distribution has been derived in Ref. [8]. Here, we use a slightly modified equation compared to their Eq. (41), since the derivation of that formula involves a multiplication with the function A (defined below). However, this function can become zero for certain complex frequencies, and multiplying the dispersion relation by A thus introduces spurious solutions. Instead of solving their Eq. (41), we directly set to zero the determinant of the matrix in their Eq. (C15), where using the same notation as in Ref. [8]. A Python program is used to solve for the complex frequencies fulfilling this dispersion relation, with very similar algorithms as in the recently introduced HYDROS code for the hybrid-kinetic system (see next section).

The hybrid-kinetic system of equations
The second model that will be examined in this work, is the hybrid kinetic-ion/fluidelectron model (which we will simply call "hybrid-kinetic" in the following). The equations of this model are obtained by taking a nonrelativistic limit of the fully kinetic equations, and taking the electron mass to zero. Retaining only a singly charged ion species, such that n i = n e , then leads to a system of equations consisting of the ion Vlasov equation, and an Ohm's law determining the electric field, which reads with n i u i = V vf i d 3 v, the resistivity η, and the electron pressure gradient ∇P e = C∇n γ e . The electromagnetic fields are constrained by Faraday's law and the pre-Maxwell version of Ampere's law (which, due to the absence of the displacement current, implies quasi-neutrality) These equations contain the full kinetic ion physics including wave-particle interactions such as Landau and transit-time damping, as well as cyclotron resonances. On the other hand, electrons appear only as a neutralizing, massless background. The equation of state for the electron pressure gradient is determined by P e = Cn γ e , with C = n 1−γ e T 0e and, in this work, γ = 1, corresponding to an isothermal equation of state.
A dispersion relation for this kind of model, valid for arbitrary propagation angle and bi-Maxwellian plasmas, has recently been derived [30], and the numerical dispersion solver HYDROS (previously employed in Ref. [31]) will be used in this work to obtain the wave solutions for that model.

The fully kinetic system of equations
Finally, the reference for this comparative analysis is provided by the nonrelativistic, fully kinetic system of equations, with one Vlasov equation (Eq. 3) for each species σ, and the full Maxwell equations including the displacement current , and the integral running over the full velocity space V. The inclusion of the displacement current in Ampére's law introduces (in normalized units) a new dimensionless parameter v A /c (with v A = B/ √ 4πm i n i ), which is set to 0.01 throughout this work, small enough to avoid discrepancies due to quasineutrality violations (in the solar wind, e.g. from Ref. [32], v A /c ≈ 2 · 10 −4 ). The fully kinetic wave solutions are obtained here using the DSHARK dispersion relation solver [33] in the limit of an isotropic Maxwellian background distribution.

Comparative study of linear wave physics
In the main part of the paper, a comparison of gyrokinetic (GK), hybrid-kinetic (HK), and fully kinetic (FK) wave physics is provided for some of the waves that are encountered in conditions found in the solar wind or magnetospheric plasmas, and which have been the center of various previous studies, namely the kinetic Alfvén wave (KAW), and the fast magnetosonic mode/whistler.

On the choice of physical parameters
For the analysis detailed in this section, the focus is on parameters that are suitable for magnetized plasmas, specifically for the small (kinetic) scales close to or below the ion gyroradius (kρ i ∼ 1, where ρ i = m i v th,i c/eB, and v th,i = 2T i /m i is the ion thermal velocity ‡.) or ion skin depth scale (kd i ∼ 1, where d i = c/ω pi ), where kinetic effects become significant. Owing to the nature of the MHD cascade [34,35], if there is a sufficiently broad spectral range between an approximately isotropic injection scale and the scale of interest, the turbulent fluctuations in said range of interest will exhibit anisotropic wavevectors with k k ⊥ , low frequency ω Ω ci , and small amplitude δB B 0 . These are the requirements for the validity of gyrokinetic theory (see Eq. 1), and both theoretical arguments and spacecraft observations [32,[36][37][38] indicate that these requirements are fulfilled in the kinetic wavenumber range of the solar wind close to 1 AU, although other observations point to the occurrence of frequencies significantly larger than Ω ci [39]. In this work, we will remain agnostic toward these questions, and will compare the three different models both in regimes that fulfill the above requirements as well as ones that do not.
Most significantly, we emphasize that all wavenumber spectra in this work are created using a constant propagation angle, which is most likely not realized in actual ‡ Note that this definition agrees with the one of Ref. [8], but differs by a factor √ 2 from the definition used in Ref. [15].
plasmas. In fact, following the critical balance arguments by Goldreich/Sridhar (GS) [34], but also Boldyrev [35], MHD turbulence is characterized by k ∝ k α ⊥ , i.e. the ratio between k and k ⊥ , and thus the propagation angle, is not independent of the wavenumber. The exponent α takes a value of 2/3 for GS, and 1/2 ≤ α ≤ 2/3 for Boldyrev, and can assume even smaller values (1/3) in the kinetic range, or for turbulence weakened by dissipation [40]. In the present work, however, our aim is not to reproduce realistic spectra, but to compare the three different models, justifying the choice of a simple, fixed, k /k ⊥ ratio.
Spacecraft observations from Refs. [32,39] indicate that the average propagation angle in the solar wind is 87.8/87.7 • , respectively, with deviations of at most ∼ 15 • [32] (∼ 30 • [39]). Therefore, most tests in the following will be performed at angles close to the average value, although a more extreme propagation angle of 60 • will be analyzed as well.
Finally, we note that the solar wind plasma usually exhibits temperatures that are anisotropic with respect to the mean magnetic field (see, e.g., Ref. [41]). Here, we focus on an isotropic background and leave a comparison of anisotropy driven instabilities (such as the mirror and the firehose mode) for future work. Previous work indicates that standard gyrokinetics should be able to reproduce the behavior of the mirror mode for almost perpendicular propagation, but may require generalizations to treat the firehose instability [42]. The GK dispersion relation used here [8] is formulated for isotropic background temperature, and would thus have to be generalized for a future comparison.

Kinetic Alfvén waves
The first test case for which the three candidate models will be compared is the kinetic Alfvén wave, i.e. the continuation of the MHD shear Alfvén wave into the range where kρ i ∼ 1 (where FLR effects become important) or kd i ∼ 1, where ion inertia becomes significant and the electron and ion motion thus decouple.

Propagation angle dependence
First, a set of fixed wave propagation angles is chosen, and scans over the magnitude of the wave vector are performed in order to obtain a dispersion relation of the KAW for this parameter set. For this analysis, β i = β e = 1 is chosen (with β σ = 8πn σ T σ /B 2 ), and wavenumbers are scanned from kd i = kρ i = 0.01 to kd i = 20.
Let us analyze the results presented in Fig. 1 first. For this figure, we set θ = 87.5 • . As is apparent from the left hand part a) of the figure, the frequencies for all models agree very well over most of the wavenumber range, 0.01 ≤ kd i 6, until the KAW reaches the cyclotron frequency. At this point, the gyrokinetic model, lacking any physics effects related to the cyclotron motion, continues to exhibit an increasing frequency, whereas the frequencies from both the hybrid-kinetic and the fully kinetic code roll over (though not exactly at the same value) and asymptote against a frequency close to the ion cyclotron frequency Ω ci .  In the right hand part, Fig. 1b), the corresponding damping rates are presented. Here, relatively good agreement between all models is observed in the wavenumber range 0.01 ≤ kd i ≤ 1. While the GK model agrees very well with full kinetics in that range, the hybrid-kinetic model underestimates the damping rates there by about 25%. This discrepancy is the first hint of a common theme that will be explored more thoroughly in the following sections, namely the absence of electron Landau and transit time damping in the hybrid-kinetic model. Importantly, and perhaps counterintuitively for a model that includes full kinetic ion physics, this can and will lead to discrepancies on ion spatial scales.
Here, this discrepancy becomes more obvious as the wavenumber is increased to kd i 1, where the gap between the HYDROS and DSHARK damping rates widens to about three orders of magnitude, until cyclotron damping becomes significant at kd i ∼ 6, closing the gap in damping rates and reestablishing agreement between the hybrid and fully kinetic physics. GK, on the other hand, missing the cyclotron damping effect, does not consistently agree with full kinetics in that range.
Before examining such discrepancies in more detail, the linear KAW dispersion relation is studied for one more propagation angle, namely θ = 60 • . The results for this scan are shown in Fig. 2. As before, very good agreement of all frequencies is observed at low wavenumbers, until the wave frequency approaches the ion cyclotron frequency at kd i ≈ 0.8. At higher k, the frequencies in the hybrid-kinetic and fully kinetic solver roll over, while in GK they pass through the cyclotron frequency and continue with an ω ∝ k 2 relation.
For this propagation angle, the low-k damping rates between all models agree very well, indicating that ion Landau damping is dominant for these parameters. At kd i ≈ 0.6,  the damping rates returned by the GK solver start to deviate from those of HYDROS and DSHARK, due to the missing cyclotron damping.
However, it is noteworthy that any agreement is found at all between GK and the other models for this propagation angle, considering that the ordering k k ⊥ is a fundamental requirement of gyrokinetics. For θ = 60 • , however, k ≈ 0.58k ⊥ , technically violating the aforementioned ordering. This finding (if found to be independent of physics parameters such as β i and β e ), is of relevance to systems like the solar wind considering the limited variability of the propagation angles found in the solar wind [32,39]. Since GK can apparently still predict useful KAW wave frequencies and damping rates for propagation angles as small as 60 • , we may expect that, if the model fails in systems like the solar wind, it is more likely to do so because of the lack of fast waves and cyclotron effects, not because of the k ordering.

Beta dependence
In this section, we compare the GK, hybrid-kinetic and fully kinetic solvers for different beta values. A fixed propagation angle of θ = 85 • is chosen (i.e. slightly less oblique than the average observed angle in the solar wind), and wavenumber spectra for two separate β i = β e values are analyzed.
First, let us examine the plots shown in Fig. 3, which were obtained for β i = β e = 5. Here, remarkably, very good agreement of wave frequencies is found for all models, across the complete wavenumber range. In this case, the KAW has a right-handed polarization and is thus not affected by the ion cyclotron resonance. The GK dispersion relation thus remains valid significantly beyond the ion cyclotron frequency. This behavior of the KAW for strongly oblique propagation has been both numerically observed and analytically derived before, and is discussed in Refs. [43,44].
Even more surprising is the comparison of the damping rate curves shown in Fig. 3b), where the GK model is found to reproduce the KAW damping rates more accurately  than the hybrid-kinetic model. At low wavenumbers kd i 0.7, all three models agree very well. At larger wavenumbers, though, electron Landau damping appears to become significant. Collisionless damping mechanisms involving the electron species are, however, completely absent from the hybrid-kinetic model, resulting in damping rates (potentially caused by anomalous ion cyclotron damping [45]) that are about one order of magnitude smaller than those of GK or full kinetics in the range kd i 1.
These circumstances need to be interpreted cautiously, however: while for these parameters the representation of the KAW is indeed more accurate in GK than in hybrid-kinetic, this may not automatically carry over to a turbulent state, where a KAW cascade may transfer energy into other waves that occur close to the cyclotron frequency, like ion Bernstein modes. The latter type of wave is not contained in gyrokinetics, and such effects are thus absent. The question about how important those effects are, cannot be answered in the linear framework of this study, and will have to be addressed in nonlinear, turbulent simulations.
As discussed in Section 3.2.1 (for a propagation angle of 87.5 • ), for β i = β e = 1 the situation is qualitatively different, as the KAW wave is left hand polarized, thus enabling the ion cyclotron resonance, and leading to results similar to the ones shown in Fig. 1. For this lower β, we also begin to observe a discrepancy of the hybrid-kinetic damping rates compared to the other two models, on ion spatial scales.
Decreasing β further to β i = β e = 0.2, we obtain the dispersion relations presented in Fig. 4. For the wave frequencies, the picture is qualitatively identical to that of Fig. 1 (β = 1, θ = 87.5 • ). The damping rate curves of Fig. 4b), however, reveal a more severe discrepancy between the hybrid-kinetic and the two other models: Across the whole ion range, 0.01 ≤ kd i 1, the hybrid-kinetic model underestimates the KAW damping rate by roughly one order of magnitude. At higher wavenumbers, where the ion Landau damping diminishes, this gap widens further, until ion cyclotron damping becomes dominant at kd i ≈ 4.  Since the discrepancy in the damping rates clearly depends on the beta parameter, another scan is performed in order to characterize this effect more stringently, this time taking a fixed wavenumber of kd i = 0.1, while varying β i = β e in the range 0.01 ≤ β i = β e ≤ 50. The wavenumber kd i is deliberately chosen to lie in the ion range of spatial scales, since this is the range where one would naively expect a model with fully kinetic ion treatment to agree with a fully kinetic model. The damping rates from this scan are shown in Fig. 5a), and their ratios from the reduced models compared to the fully kinetic result are plotted in Fig. 5b). For these parameters, GK agrees very well with full kinetics across the whole beta range. On the other hand, consistent with the picture that emerged above, the hybrid-kinetic model shows excellent agreement at high β 2, while for lower beta values a discrepancy arises. This discrepancy is about 25% for β = 1, about one order of magnitude for β = 0.2, and three orders of magnitude for β = 0.1, and increases quickly for lower β.
A physical interpretation of this discrepancy may be obtained by considering that for an Alfvén wave in this range The ion Landau resonance occurs for ions traveling at speeds similar to the propagation velocity of the Alfvén wave, i.e. v ≈ v A . For a Maxwellian distributed plasma, such particles are most abundant for β i values such that v A v th,i , i.e. the ion Landau resonance is most effective for β i 1, in agreement with Fig. 5. On the other hand, for β i 1 the Alfvén wave travels at a velocity faster than most ions and thus detunes from the ion Landau resonance, resulting in the strongly reduced damping rates observed in the hybrid model.
In the case of electron Landau damping, the above statements apply in a similar way, except that v A v th,e is now the relevant condition. Since in a hydrogen plasma the electron thermal velocity is roughly a factor 43 faster than the ions', this condition is fulfilled for a much broader range of β values, making the electron Landau damping more resilient against variations of this parameter. For the examined β values, it is thus the electron Landau damping that remains once the ion Landau damping is removed and that explains the discrepancy between the hybrid-kinetic and the fully kinetic model.

Mass ratio effect
Motivated by the above results, we examine in this section the consequences of a reduced mass ratio on the behavior of a kinetic Alfvén wave. This question is of interest, as fully kinetic PIC or Eulerian turbulence simulations are often extremely challenging in a computational sense, enforcing either a reduced dimensionality of the simulations, or the removal of physical scale separations such as those between ions and electrons. The latter is achieved by reducing the mass ratio between the species, leading to a compression of their spatiotemporal scales. In order to examine the effect of such a reduced mass ratio on the kinetic Alfvén wave, we choose once more a propagation angle of θ = 85 • , and we perform a wavenumber scan for β i = β e = 1. The results of this scan are depicted in Fig. 6. As can be observed in Fig. 6a), the wave frequencies are not significantly altered by changing the mass ratio, although small differences are visible close to the ion cyclotron frequency. On the other hand, the damping rates plotted in Fig. 6b) show a visible discrepancy (about 65% overprediction for reduced mass ratio) across the ion range of wavenumbers, and more significantly so in the range 1 kd i 4, where ion Landau damping is weakened and ion cyclotron damping is not yet active.
Considering the findings of Sec. 3.2.2, we must suspect a beta dependence of this finding, and based on the discussion at the end of that section, we may expect to find more significant differences at lower β, where electron Landau damping dominates. In  Another finding is related to the switching of the KAW polarization that was observed in Fig. 3 for β = 5. The exact point in parameter space where this switching occurs depends on the ion/electron mass ratio, so that different results are obtained both for frequencies and damping rates in that regime, depending on the mass ratio employed. In Fig. 7b), another wavenumber scan is shown to illustrate this for β = 5, where different mass ratios yield damping rate differences up to a factor 4.5 at large kd i . These discrepancies occur because the wavenumber where ion cyclotron damping becomes dominant shifts with mass ratio, exposing the variation in (the otherwise subdominant) electron Landau damping.

Fast magnetosonic waves/Whistler modes
In this section, the focus lies on the fast magnetosonic mode and its potential transition to a whistler mode which exhibits an ω ∝ k 2 scaling. This wave is ordered out of gyrokinetics, reducing the set of models to just the hybrid-kinetic and the fully kinetic theory. The procedure adopted here is similar to that for the kinetic Alfvén wave, and we first analyze the dependence of the dispersion relation on the wavenumber, for two different propagation angles. relation until the ion cyclotron frequency is reached. For the chosen propagation angle, this wave experiences ion cyclotron damping and stays below the cyclotron frequency. Despite the qualitative agreement, we note that a shift between the two frequency curves can be observed in the range 0.01 ≤ kd i 0.6, until the ion cyclotron frequency is reached. This deviation was found before in Ref. [30] and could be traced back to the massless electron approximation underlying the hybrid-kinetic model. Conversely, it is possible to obtain the same result from the fully kinetic dispersion solver by numerically approximating the massless electron limit.
In the same wavenumber range, there is a striking difference in the measured damping rates: while DSHARK reports a rather strongly damped fast mode (|γ| /ω 0.1), HYDROS finds a completely undamped mode (within the numerical accuracy of the code). With increasing wavenumber, agreement is only obtained when the ion cyclotron frequency is approached, and for kd i 1 the codes then agree very well on both frequencies and damping rates.
This severe discrepancy at low wavenumbers (once more, on ion scales) is due to electron transit time damping (i.e. caused by the near-constant mirror force observed by an electron that travels at similar speed to the fast wave) [46]. This effect is not contained in the hybrid-kinetic model solved by HYDROS, explaining the observed deviations. In contrast to the kinetic Alfvén wave, this effect plays an important role even for β = 1, which can be explained by the larger phase velocity of the fast mode compared to the KAW: in order to obtain a resonant behavior with electrons, the wave needs to travel at a velocity v ph v th,e . The phase velocity of the fast wave for these parameters is roughly given by ω ≈ √ 3k ⊥ v A ≈ 40k v A = 40k v th,i ≈ k v th,e (using Eq. (14) from Ref. [43]), in good agreement with the above condition. Interestingly, when the wave reaches the cyclotron frequency, its linear frequency dependence on k is lost, leading to a detuning of the resonance and a relatively sharp reduction of the damping rate, above kd i ≈ 0.6. Next, let us study the fast magnetosonic mode/whistler dispersion for a propagation angle of 60 • , keeping β i = β e = 1. The results of these scans are presented in Fig. 9. For these parameters, very good quantitative agreement of the wave frequencies between the hybrid-kinetic and fully kinetic solvers is found. In addition, the fast wave now passes smoothly through the ion cyclotron resonance, and for wavenumbers close to kd i ≈ 3, the dispersion relation transitions from its linear (in k) to the well-known quadratic whistler behavior. Although not shown in more detail here, this behavior of the wave is found (for β = 1) for propagation angles θ 63 • . For more oblique angles, the fast magnetosonic wave experiences ion cyclotron damping and ceases to propagate at the ion cyclotron frequency. Note that for such parameters high-k, high-frequency waves with an ω ∝ k 2 dependence are still found, but they do not smoothly connect to the sub-Ω ci fast magnetosonic mode.
Examining in turn the linear damping rates, again a region is found where electron transit time damping dominates (kd i 0.4), leading to a two order of magnitude discrepancy between HYDROS and DSHARK in that regime. At higher wavenumbers, even though the wave does not have a cut-off at Ω ci , it still experiences enhanced damping rates between 0.4 kd i 3. This damping is ion cyclotron damping and is thus captured well by HYDROS. At still higher wavenumbers, electron damping dominates again, widening the gap between the hybrid-kinetic and fully kinetic damping rates again. Finally, we note that, although not shown here, for even less oblique propagation (i.e. moving further away from observed solar wind propagation angles) ion Landau damping becomes increasingly important, thus improving the agreement between the two models significantly.
3.3.2. Beta dependence Following the strong discrepancy between completely undamped (hybrid-kinetic) and rather heavily damped (fully kinetic) fast modes in the previous section, we will now characterize the parameter space for which this finding is of relevance. Indeed, motivated by the findings regarding the beta dependence of electron Landau damping of kinetic Alfvén waves, the same kind of analysis is now performed for the fast magnetosonic mode, as similar arguments may be expected to hold for this kind of wave.
As before in Sec. 3.2.2, the propagation angle is fixed to θ = 85 • and wavenumber scans are performed for various values of β i = β e . In Fig. 10, only damping rates are shown, since the frequencies generally exhibit satisfactory agreement between the hybrid-kinetic and fully kinetic models. In Fig. 10a), a high beta example for β = 20 is chosen, remembering that for KAWs ion Landau damping was found to dominate in that parameter regime, masking the lack of electron physics in the hybrid model. For the fast magnetosonic mode, however, even for the quite high β = 20 there is a significant wavenumber region up to kd i ≈ 0.15 where the wave is undamped in the hybrid-kinetic model, in contrast to the fairly strong damping (|γ| /ω ≈ 0.1) that occurs in the fully kinetic system. At larger wavenumber, cyclotron damping dominates, which is in turn well matched by the hybrid-kinetic model. This picture remains qualitatively the same for a wide range of β values, although the curves shift (in d i normalization) in wavenumber space. Fig. 10b) depicts a similar wavenumber scan for a low β = 0.01. Here, because of the lower wave frequency of the fast magnetosonic mode at low kd i a very wide undamped wavenumber range is found in the hybrid-kinetic model, with cyclotron damping occurring only at wavenumbers of kd i 10. The mode is more weakly damped at this low β value, but the qualitative picture from before remains intact.
The very wide β range across which electron transit time damping is found to dominate the picture here can be explained, as was hinted before, by the much wider electron velocity distribution, owing to their large thermal velocity. Heuristically speaking, when varying β, a wave is detuned much more easily from an ion (Landau or transit-time) resonance than from an electron resonance.

Mass ratio effect
Finally, we now analyze the effect of a reduced mass ratio on the fast magnetosonic mode, using the DSHARK solver. The findings of the previous section, where a massless electron description was found to result in practically undamped fast modes, suggest that a mass ratio reduction may have a similarly strong impact. Here, once more wavenumber scans for the mass ratios m i /m e = 1836, 400, and 100 are performed, for β i = β e = 1 and a propagation angle of 85 • .
For these parameters, similar to the findings for the KAW, there is no strong impact of the mass ratio on the wave frequency, see Fig. 11a). However, on ion scales below kd i ≈ 0.6 the wave damping rate (Fig. 11b)) turns out to be sensitive to the mass ratio: while for a mass ratio of m i /m e = 400 the damping rates deviate only be a few percent, the simulations with a mass ratio of 100 underestimate the real damping rates by a full order of magnitude. This behavior can again be understood in terms of the detuning of the resonance as the electrons are made heavier and heavier -for a mass ratio of 100, the fast magnetosonic mode can interact only with the weakly populated tail of the heavy electron distribution. At wavenumbers of kd i 0.6 on the other hand, ion cyclotron damping starts to dominate, which is not affected by the reduced mass ratio.
In Fig. 12a), the ratios of the damping rates obtained for reduced mass ratio to the real damping rates are plotted. As can be seen, a mass ratio of 400 gives a reasonable approximation to the real damping rates, except in the region where the transition between the two different damping mechanisms occurs. In contrast to that, for a mass ratio of 100 the real damping rate is severely underpredicted for low wavenumber and is only recovered with reasonable accuracy in the range of dominant ion cyclotron damping, at large wavenumber.
Finally, in Fig. 12b), these damping rate ratios are plotted for fixed wavenumber kd i = 0.1, while varying beta over several orders of magnitude. In simulations with a  mass ratio of m i /m e = 400, the fast wave damping rates are relatively well represented for large beta (with a discrepancy of about 30% for β 1), but reach a discrepancy of more than an order of magnitude when β < 0.1. The result for m i /m e = 100 is more concerning, however, as over the entire range of examined beta values the damping rates are underestimated by at least a factor 2.5, and by more than an order of magnitude for β < 1. The only exception occurs at β = 100, where for the fixed wavenumber of kd i = 0.1 the fast wave reaches a frequency close to Ω ci , so that ion cyclotron damping dominates here, which does not depend on the mass ratio.

Summary and conclusions
In the present study, the linear wave physics contained in the gyrokinetic, the hybridkinetic (combining a fully kinetic ion treatment and fluid electron model) and the full kinetic model of plasma physics were analyzed and compared. For this purpose, we focused on two different wave modes, namely the kinetic Alfvén wave and the fast magnetosonic/whistler wave. For both kinds of wave the dispersion relations were studied for two different propagation angles, one of them chosen to be 87.5 • , close to the average propagation angle found in solar wind plasmas, and the other close to the maximal observed deviation from the aforementioned average angle, namely to 60 • [39].
For the kinetic Alfvén wave (KAW), it was found that the gyrokinetic model (GK) generally agrees very well with full kinetics as long as the ion cyclotron frequency Ω ci is not reached. In some cases, this is even true for frequencies much higher than Ω ci , provided that the KAW is right-hand polarized. However, energy transfer between KAWs and other waves present at such high frequencies (that are missing from GK) cannot be accounted for within the gyrokinetic model.
As expected, the hybrid-kinetic model was found to agree very well with the fully kinetic model, as long as electron wave-particle interactions do not dominate the wave damping. If they do, however, the hybrid-kinetic model often severely underestimates the linear wave damping rates, and, perhaps counterintuitively so, even on ion spatial scales. Since the phase velocity of the KAW is close to the Alfvén velocity (and v A ∝ v th,i / √ β i ), the KAW detunes from the ion Landau resonance at low β. Then, electron Landau damping is the dominant damping mechanism, resulting in an underprediction of damping rates by the hybrid-kinetic model.
The fast magnetosonic wave, on the other hand, is ordered out from GK (even in cases where its frequency is below Ω ci ), so only the hybrid-kinetic and fully kinetic models could be compared here. For a propagation angle of 87.5 • and β i = β e =1, a striking disagreement between the hybrid-kinetic and the fully kinetic model is found: while the wave is rather strongly damped in the latter case, the hybrid model finds a completely undamped mode. In this case, the discrepancy is caused by the lack of electron transit time damping in the hybrid model. Even for less oblique propagation (θ = 60 • ), the damping rates from both models still disagree by almost two orders of magnitude on ion scales, although ion Landau damping then starts to become more important. The beta dependence of this effect was studied, showing that this observation is robust across a wide range of beta values, owing to the broad (in velocity space) electron transit time resonance.
Motivated by the fact that many fully kinetic simulations are carried out with reduced mass ratio in order to make these simulations feasible, the fully kinetic model with reduced mass ratio was benchmarked against its real-mass ratio counterpart. Using a mass ratio of m i /m e = 100, the effect on KAWs was relatively moderate and mostly limited to β 1. The fast magnetosonic mode, on the other hand, was more gravely affected and it was found that ion scale damping rates of these modes are underestimated by at least a factor 2.5 across the whole studied beta range, and up to several orders of magnitude for β 1. While 3D fully kinetic simulations using the real proton/electron mass ratio will remain very demanding for the time being, the findings of this work should be of use for interpreting existing and future simulations, and for obtaining projections to the real systems they are meant to describe.
Finally, we would like to remark that the studies performed in this work have barely scratched the surface of the comparisons that are possible. Given the popularity of both gyrokinetic and hybrid-kinetic simulations, we believe that the availability of easy-to-use dispersion solvers for both these models is essential for a realistic assessment of the models and the simulations performed with them. For this reason the hybrid-kinetic solver HYDROS has been made available on Github [47], as was the DSHARK solver before [48].