1 Introduction
One of the distinctive characteristics of plasmas is that they have time and length scales that are typically well separated. For example, in an electron–proton plasma the mass ratio between the two species is approximately equal to 1836, and this determines a separation of the same factor in the respective cyclotron frequencies and of a factor
${\sim}\sqrt{1836}\simeq 43$
in the respective Larmor radii and plasma frequencies. Such a large separation of scales translates very often into unaffordable costs, from a computational point of view. In particular, it is well known that a fully kinetic treatment, where two Vlasov equations (one for each species) are coupled to Maxwell’s equations, is prohibitively expensive in multiple dimensions. For this reason, several reduced models have been studied in the past that are able to retain some of the kinetic features of the plasma at a reduced computational cost (Matthews Reference Matthews1994; Waltz et al.
Reference Waltz, Staebler, Dorland, Hammett, Kotschenreuther and Konings1997; Cheng & Johnson Reference Cheng and Johnson1999; Park et al.
Reference Park, Belova, Fu, Tang, Strauss and Sugiyama1999; Goswami, Passot & Sulem Reference Goswami, Passot and Sulem2005; Brizard & Hahm Reference Brizard and Hahm2007; Valentini et al.
Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007). Such models could be ordered in a hierarchy that transitions all the way from kinetic to fluid theory (Dendy Reference Dendy1995; Tronci et al.
Reference Tronci, Tassi, Camporeale and Morrison2014; Tronci & Camporeale Reference Tronci and Camporeale2015; Camporeale et al.
Reference Camporeale, Delzanno, Bergen and Moulton2016). In this paper, we focus specifically on two models that have been extensively used for solar wind turbulence studies: the hybrid (ion kinetic and electron fluid) and the gyrokinetic models. The philosophy behind these two models is quite different. In the hybrid model the heavy species (protons) is treated kinetically and its distribution function obeys the Vlasov equation (which is in principle six-dimensional), while the light species (electrons) is treated as an inertialess neutralizing fluid. Thus, the computational saving is achieved by factoring out from the model the time and spatial scales associated with electrons. We note that the term ‘hybrid’ is used widely in different contexts with different meanings, but here we use it with the specific meaning given. On the other hand, the gyrokinetic model treats both species kinetically, but it performs a time average on the ion-cyclotron period. As a result, the six-dimensional distribution functions collapse into five dimensions. Interestingly, such dimensionality reduction of the particle distribution functions (and the fact that the dynamics faster than the ion gyroperiod is neglected) results in a sizable saving in computational time.
Like all reduced models, the derivation of the hybrid and gyrokinetic models is based on a certain number of assumptions. They both make use of the so-called quasi-neutrality condition, which asserts that ion and electron densities are equal,
$n_{i}=n_{e}$
. This is satisfied as long as the Debye length is much smaller than the characteristic spatial scales, and the characteristic frequencies are much smaller than the electron plasma frequency. As a consequence, the displacement current in Ampere’s law is neglected in both models. Furthermore, the gyrokinetic theory assumes that frequencies are much smaller than the ion-cyclotron frequency,
$\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D6FA}_{i}$
, and that the wave spectrum is very anisotropic,
$k_{||}\ll k_{\bot }$
. The hybrid model, instead, necessitates of a closure in terms of an equation of state for the electrons.
The fundamental problem that we tackle in this paper is: what does much smaller mean?
Once a model (and the computer code that solves it) has been extensively used and benchmarked in a certain regime, it is certainly tempting to push the simulations into unexplored regimes. Although there are cases for which the reduced models are obviously inadequate, there exist a large range of ‘grey areas’, where their validity is not obvious. In a sense, the issue of validating a reduced model is an ill-posed problem. In principle, the validation of the results in these grey areas can only be performed by comparing with a fully kinetic model. On the other hand, if such solutions from the fully kinetic models were available, one would not need to use a reduced model!
The rationale for this work stems from the fact that recent simulations performed with hybrid and the gyrokinetic codes might be considered the current state of the art in studies concerned with solar wind turbulence (see, e.g. Valentini, Califano & Veltri Reference Valentini, Califano and Veltri2010; Howes et al. Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011; Servidio et al. Reference Servidio, Valentini, Califano and Veltri2012; Franci et al. Reference Franci, Verdini, Matteini, Landi and Hellinger2015). Indeed, it is now widely recognized that a kinetic treatment is necessary for studying solar wind turbulence at small scales (Camporeale & Burgess Reference Camporeale and Burgess2011; Sahraoui, Belmont & Goldstein Reference Sahraoui, Belmont and Goldstein2012; Boldyrev et al. Reference Boldyrev, Horaites, Xia and Perez2013; Haynes, Burgess & Camporeale Reference Haynes, Burgess and Camporeale2014; Vásconez et al. Reference Vásconez, Valentini, Camporeale and Veltri2014; Haynes et al. Reference Haynes, Burgess, Camporeale and Sundberg2015). Therefore, it is legitimate to ask what is the range of validity of such models and when do the results start to deviate considerably from the correct solutions. Here, we focus on comparing solutions obtained within the linear theory. Clearly, the ability of a model to describe the property of linear waves, within a certain accuracy, can be interpreted as a necessary (but not sufficient) condition for the model to yield reliable results in the nonlinear regime.
The scope of this paper is twofold. First, we directly compare linear theory solutions obtained from the hybrid and the gyrokinetic models with the complete Vlasov–Maxwell model, for wavevectors
$k$
ranging from
$k\unicode[STIX]{x1D70C}_{i}=0.1$
to
$k\unicode[STIX]{x1D70C}_{i}=44$
(with
$\unicode[STIX]{x1D70C}_{i}$
the ion Larmor radius), and for a wide range of propagation angles
$\unicode[STIX]{x1D703}$
, ranging from
$10^{\circ }$
to
$85^{\circ }$
; second, we present a robust and novel methodology to perform such a comparison in a fair way, that takes into account the fact that for large wavevector and/or propagation angles, obtaining linear solutions is a non-trivial numerical exercise due to the presence of many normal modes.
The paper is organized as follows. We introduce the methodology and the parameters chosen for a test case in § 2, where we also define the errors that will be measured. Section 3 shows the comparison between hybrid, gyrokinetic, and fully kinetic solutions for a wide range of cases. Conclusions are drawn in § 4.
2 Methodology and chosen parameters
We study a plasma composed of electrons and protons, with the mass ratio set to the physical value
$m_{i}/m_{e}=1836$
, in Cartesian coordinates
$(\hat{i} ,{\hat{j}},\hat{k})$
. As is customary in linear theory, we study a uniform Maxwellian plasma with density
$n_{0}$
in a constant magnetic field
$\boldsymbol{B}=B_{0}\hat{k}$
, that is perturbed with a small-amplitude plane wave of the form
$\exp [\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D714}t)]$
, with wavevector
$\boldsymbol{k}=(k_{\bot },0,k_{||})$
and complex frequency
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+i\unicode[STIX]{x1D6FE}$
.
We focus on a single choice of parameters that are relevant for solar wind turbulence studies. In order to facilitate comparisons with published literature, we use the typical solar wind parameters studied in Salem et al. (Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012):
$B_{0}=11$
nT,
$T_{e}=13.0$
eV (electron temperature),
$T_{p}=13.6~\text{eV}$
(proton temperature),
$n_{0}=9~\text{cm}^{-3}$
. The proton plasma beta is
$\unicode[STIX]{x1D6FD}_{p}=0.4$
and the electron plasma to cyclotron frequency ratio is
$\unicode[STIX]{x1D714}_{pe}/\unicode[STIX]{x1D6FA}_{e}\simeq 88$
. The ratio between Alfvén and ion thermal velocities is
$v_{A}/v_{ti}=1.57$
, and the ratio between the Alfvén velocity and the speed of light is
$v_{A}/c=2.6\times 10^{-4}$
.
In all of the three cases considered (fully kinetic, hybrid and gyrokinetic), the linear wave solutions are obtained by finding the complex frequencies
$\unicode[STIX]{x1D714}$
for which the determinant of a
$3\times 3$
matrix
$\unicode[STIX]{x1D63F}$
vanishes.
$f(\unicode[STIX]{x1D714},\boldsymbol{k})=\det (\unicode[STIX]{x1D63F})$
is a nonlinear function of both
$\boldsymbol{k}$
and
$\unicode[STIX]{x1D714}$
, and the standard way to numerically solve the problem, for a fixed
$\boldsymbol{k}$
, is by employing an iterative root finder based on Newton’s method. Unfortunately, although it is straightforward to extend Newton’s method for finding zeros in the complex plane, the method can be very sensitive to the choice of an initial guess. Moreover, due to the possible presence of multiple solutions which are close to each other in the complex plane, it is not always easy to track a solution
$\unicode[STIX]{x1D714}$
varying
$\boldsymbol{k}$
. In other words, the solver can easily jump from one solution branch to another whenever two solutions become close enough. The paper by Podesta (Reference Podesta2012) explains how to reduce the numerical error in calculating the elements of
$\unicode[STIX]{x1D63F}$
in the case of the Vlasov–Maxwell system. Podesta (Reference Podesta2012) also suggests identifying possible roots by visual inspection of the function
$1/|f|$
, where the peaks might indicate roots. An alternative way, which we think is more robust, less time consuming and can be easily automated, is to analyse separately the real
$(\text{Re})$
and imaginary
$(\text{Im})$
parts of the function
$f(\unicode[STIX]{x1D714})$
, and to identify the curves in the complex plane, where
$\text{Re}(f)=0$
and
$\text{Im}(f)=0$
, that is where they change sign. Clearly, the intersections of such curves (
$\text{Re}(f)=\text{Im}(f)=0$
) are the roots
$f(\unicode[STIX]{x1D714})=0$
that can be used as initial guess for the Newton solver to achieve a better accuracy. We have employed this method to identify all the roots in a given region of the complex plane, in § 2.1.1
We refer to Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006) for a derivation of the gyrokinetic linear theory. The linear theory for the hybrid model can easily be derived by noting that, in the standard derivation of the Vlasov–Maxwell linear theory, the dielectric tensor
$\unicode[STIX]{x1D73A}(\unicode[STIX]{x1D714},\boldsymbol{k})$
is additive by species (see, e.g. Stix Reference Stix1962):
$$\begin{eqnarray}\unicode[STIX]{x1D73A}(\unicode[STIX]{x1D714},\boldsymbol{k})=\mathbf{1}+\mathop{\sum }_{s}\unicode[STIX]{x1D712}_{s}(\unicode[STIX]{x1D714},\boldsymbol{k}),\end{eqnarray}$$
where
$\mathbf{1}$
is the unit dyadic,
$\unicode[STIX]{x1D712}_{s}$
is the susceptibility of species
$s$
and the sum is over plasma species. Therefore, in a hybrid model it is sufficient to calculate the electron susceptibilities
$\unicode[STIX]{x1D712}_{e}$
derived within a two-fluid linear model. Here, we have employed the so-called warm plasma dispersion relation derived in Swanson (Reference Swanson2012), which assumes the equation of state
$p_{e}/n_{e}=v_{te}^{2}$
(with
$p_{e}$
,
$n_{e}$
and
$v_{te}$
the electron pressure, density and thermal velocity, respectively). Moreover, in order to enforce the quasi-neutrality condition, it is sufficient to suppress the unit dyadic term in the dielectric tensor, which is thus defined simply as the sum of the ion susceptibility (derived in kinetic theory) and the electron one (derived in two-fluid theory).
2.1 Diagnostics
The results obtained from hybrid and gyrokinetic linear theory are compared with the ones obtained in the fully kinetic Vlasov–Maxwell theory. We divide the parameter space into oblique propagation, for
$\unicode[STIX]{x1D703}\leqslant 60^{\circ }$
, and quasi-perpendicular propagation, for
$\unicode[STIX]{x1D703}>60^{\circ }$
. While the diagnostics used for oblique propagation are based on the comparison of single normal modes, we argue that such diagnostics are not sufficiently informative for the quasi-perpendicular regime. Typically, among the infinite number of modes which can in principle exist as solutions of the full linear Vlasov–Maxwell system, one is interested in the few that have lower damping rates, while the others (usually referred to as ‘heavily damped’) are of little interest because, supposedly, they do not play a significant role for mode coupling. We note that this view, although probably correct when the modes are well separated, is questionable when there are a multitude of normal modes clustered together. In this case, one would like to take into consideration the overall spectral properties of the linear operator. Such an approach is justified by the non-normality of the Vlasov–Maxwell linear operator. We refer to, e.g. Camporeale, Burgess & Passot (Reference Camporeale, Burgess and Passot2009), Camporeale, Passot & Burgess (Reference Camporeale, Passot and Burgess2010), Podesta (Reference Podesta2010), Friedman & Carter (Reference Friedman and Carter2014), Ratushnaya & Samtaney (Reference Ratushnaya and Samtaney2014) and citations therein, for a discussion on non-normal linear operators in plasma physics. In short, and for the purpose of this paper, we recall that non-normality is due to the non-orthogonality of the eigenvectors of a linear operator and is connected to the phenomena of transient growth, by-pass transitions and generally to mode coupling. An important characteristic of non-normal operators is that their eigenvalues are highly sensitive to perturbations in the sense that the spectrum can be easily distorted in the complex plane when the original operator is slightly perturbed. This concept gives rise to the idea of pseudospectrum, which is extensively discussed in Trefethen & Embree (Reference Trefethen and Embree2005) and, for the case of the Vlasov–Maxwell system, in Camporeale (Reference Camporeale2012). For instance, figure 1 shows a subset of the normal modes (obtained in the fully kinetic model), in the complex plane, for
$\unicode[STIX]{x1D703}=80^{\circ }$
. The damping rate
$\unicode[STIX]{x1D6FE}$
and the real frequency
$\unicode[STIX]{x1D714}_{r}$
are on horizontal and vertical axis, respectively (here and in all figures these parameters are normalized to the ion-cyclotron frequency
$\unicode[STIX]{x1D6FA}_{i}$
). The different symbols are for
$k\unicode[STIX]{x1D70C}_{i}=0.1$
(blue dots),
$k\unicode[STIX]{x1D70C}_{i}=1$
(red circles) and
$k\unicode[STIX]{x1D70C}_{i}=5$
(black crosses). Clearly, a distinction between modes purely based on their damping rate becomes less and less possible with increasing wavevector. Moreover, there is no reason to attribute a particular importance to some well-characterized modes (that is, well characterized in the low
$k$
and oblique propagation regime), such as the kinetic Alfvén wave (KAW) mode. This is particularly clear by inspection of figure 2. Here, we plot in (a) the subset of normal modes that are confined in the complex plane within the domain
$-2.5<\unicode[STIX]{x1D6FE}<0$
and
$0<\unicode[STIX]{x1D714}_{r}<20$
. The blue dots and red circles denote the fully kinetic and the hybrid solutions, respectively. The kinetic Alfvén wave, that is the root of the dispersion relation that when tracked to low
$k$
and low
$\unicode[STIX]{x1D703}$
is connected to the Alfvén wave mode, is denoted with a blue cross (
$\unicode[STIX]{x1D714}_{r}=0.47$
,
$\unicode[STIX]{x1D6FE}=-1.29$
). It is evident that, in this regime, the KAW is surrounded by many other modes, some of them with similar frequency but lower damping rate. Nevertheless, the KAW has been attributed a predominant role in some of the recent solar wind literature (Sahraoui et al.
Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009; Salem et al.
Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012). Even if one would like to adhere to the classical view where the least damped modes are more important than the heavily damped ones, it is clear that the KAW is not special. Figure 2(b) shows the same modes, but now sorted as a function of their damping rate per wave period
$(|\unicode[STIX]{x1D6FE}|/\unicode[STIX]{x1D714}_{r})$
(horizontal axis) and their polarization
$P$
(vertical axis). The polarization is here defined following Hunana et al. (Reference Hunana, Goldstein, Passot, Sulem, Laveder and Zank2013), as
$P=\arg (E_{y}/E_{x})/\unicode[STIX]{x03C0}$
, with
$E_{y}$
and
$E_{x}$
complex amplitudes of the electric field. Hence, negative/positive values of
$P$
represent left-/right-hand polarization. Interestingly, in this regime, the KAW is the most rapidly damped mode per wave period. A final remark on figure 2 is that one can see that, despite being in the electron range (i.e.
$k\unicode[STIX]{x1D70C}_{i}>1$
), the hybrid model captures most of the normal modes quite accurately. It creates a spurious mode with
$\unicode[STIX]{x1D714}_{r}=14$
and
$\unicode[STIX]{x1D6FE}=-0.08$
due to the inability to account for electron Landau damping.
Normal modes derived with the fully kinetic model. Damping rate (
$\unicode[STIX]{x1D6FE}$
) and real frequency
$(\unicode[STIX]{x1D714}_{r})$
are on the horizontal and vertical axis, respectively. The angle of propagation is
$\unicode[STIX]{x1D703}=80^{\circ }$
. Blue dots, red circles and black square denote
$k\unicode[STIX]{x1D70C}_{i}=0.1$
, 1, 5, respectively.
Normal modes for
$\unicode[STIX]{x1D703}=80^{\circ }$
and
$k\unicode[STIX]{x1D70C}_{i}=5$
for fully kinetic (blue dots) and hybrid (red circles) models. (a) Solutions in the complex plane
$(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D714}_{r})$
. (b) Solutions plotted as functions of
$|\unicode[STIX]{x1D6FE}|/\unicode[STIX]{x1D714}_{r}$
(horizontal axis) and polarization
$P$
(vertical axis). In both panels the kinetic Alfvén wave is denoted with a blue cross.
2.1.1 Definition of errors
We define the following two measures of error that assess the deviation of the frequency
$\unicode[STIX]{x1D714}$
(obtained either with the hybrid or gyrokinetic model), with respect to the correct value obtained with the Vlasov–Maxwell theory, denoted as
$\unicode[STIX]{x1D714}_{VM}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}=\frac{|\unicode[STIX]{x1D714}_{VM}-\unicode[STIX]{x1D714}|}{|\unicode[STIX]{x1D714}_{VM}|} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D73A}_{wave}=\Vert \text{Re}(\exp [-\text{i}\unicode[STIX]{x1D714}_{VM}t])-\text{Re}(\exp [-\text{i}\unicode[STIX]{x1D714}t])\Vert _{2}. & \displaystyle\end{eqnarray}$$
The error
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
in (2.2) is simply the relative distance between the two frequencies, in the complex plane, and it does not distinguish the measure in which the real and the imaginary part (damping rate) of the frequency contribute to the error. In order to calculate the error denoted as
$\unicode[STIX]{x1D73A}_{wave}$
in (2.3), we define a time interval,
$t=[0,2\unicode[STIX]{x03C0}/\text{Re}(\unicode[STIX]{x1D714}_{VM})]$
, that covers one wave period (with respect to
$\unicode[STIX]{x1D714}_{VM}$
), is discretized in 200 points and we calculate the
$L_{2}$
norm of the difference of the real part of the wave amplitude (normalized by the number of points). This metric is relevant because it comprises the information on how much a wave with frequency
$\unicode[STIX]{x1D714}$
will differ, during one wave period, from the reference wave with frequency
$\unicode[STIX]{x1D714}_{VM}$
, by modelling its actual time evolution.
Pseudospectra in the complex plane
$(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D714}_{r})$
for
$k\unicode[STIX]{x1D70C}_{i}=1$
and
$\unicode[STIX]{x1D703}=80^{\circ }$
. (a) Fully kinetic; (b) Hybrid; (c) Gyrokinetic. The colour plot represents
$f(\unicode[STIX]{x1D714})=\det (\unicode[STIX]{x1D63F})$
in logarithmic scale.
Pseudospectra in the complex plane
$(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D714}_{r})$
for
$k\unicode[STIX]{x1D70C}_{i}=1$
and
$\unicode[STIX]{x1D703}=85^{\circ }$
. (a) Fully kinetic; (b) Hybrid; (c) Gyrokinetic. The colour plot represents
$f(\unicode[STIX]{x1D714})=\det (\unicode[STIX]{x1D63F})$
in logarithmic scale.
Range of validity of the gyrokinetic model. Square and diamond symbols indicate an error in the real frequency of 10 % and 30 %, respectively. The perpendicular wavevector
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
and the angle of propagation
$\unicode[STIX]{x1D703}$
are on the horizontal and vertical axis. The model is valid above the symbols.
Range of validity of the gyrokinetic model. Dispersion relation for Alfvén/ion-cyclotron wave, plotted as
$\unicode[STIX]{x1D714}_{r}/k_{||}v_{A}$
as a function of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. The blue line denotes the gyrokinetic solution, which is independent of the angle of propagation. The red, yellow and purple lines denote the fully kinetic solution for angles of propagation
$\unicode[STIX]{x1D703}=80^{\circ },85^{\circ },89^{\circ }$
.
Dispersion relation for Alfvén/ion-cyclotron waves for angles of propagation
$\unicode[STIX]{x1D703}=10^{\circ },30^{\circ },60^{\circ }$
(curves from left to right). Blue, red and green lines are for fully kinetic, hybrid and gyrokinetic models, respectively. (a) Real frequency
$\unicode[STIX]{x1D714}_{r}$
; (b) Damping rate
$\unicode[STIX]{x1D6FE}$
.
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis. All axes are in logarithmic scale.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=10^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=30^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=60^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Dispersion relation for Alfvén/ion-cyclotron waves for angles of propagation
$\unicode[STIX]{x1D703}=80^{\circ },85^{\circ },89^{\circ }$
(curves from left to right). Blue, red, and green lines are for fully kinetic, hybrid and gyrokinetic models, respectively. (a) Real frequency
$\unicode[STIX]{x1D714}_{r}$
; (b) Damping rate
$\unicode[STIX]{x1D6FE}$
.
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis. All axes are in logarithmic scale.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=80^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=85^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for Alfvén/ion-cyclotron wave with
$\unicode[STIX]{x1D703}=89^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Dispersion relation for ion-acoustic waves for angles of propagation
$\unicode[STIX]{x1D703}=10^{\circ },30^{\circ },60^{\circ }$
(curves from left to right). Blue, red and green lines are for fully kinetic, hybrid, and gyrokinetic models, respectively. (a) Real frequency
$\unicode[STIX]{x1D714}_{r}$
; (b) Damping rate
$\unicode[STIX]{x1D6FE}$
.
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis. All axes are in logarithmic scale.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=10^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Moreover, we will use one more metric to compare results in the quasi-perpendicular regime. Here, we seek to have a measure that gives us information on the global difference between the spectra calculated with the three methods. For this purpose, we rely on the (generalized) definition of the
$\unicode[STIX]{x1D73A}$
-pseudospectrum
$\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D73A}}(\unicode[STIX]{x1D63F})$
, for a fixed
$\boldsymbol{k}$
:
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D73A}}(\unicode[STIX]{x1D63F})=\{z\in \mathbb{C}:|\text{det}(\unicode[STIX]{x1D63F}(z))|\leqslant \unicode[STIX]{x1D73A}\}.\end{eqnarray}$$
Clearly, the pseudospectrum reduces to the standard spectrum (set of eigenvalues) for
$\unicode[STIX]{x1D73A}\rightarrow 0$
. For instance, figures 3 and 4 compare the fully kinetic (a), hybrid (b) and gyrokinetic (c) pseudospectra for the case
$k\unicode[STIX]{x1D70C}_{i}=1$
and
$\unicode[STIX]{x1D703}=80^{\circ }$
(figure 3) or
$\unicode[STIX]{x1D703}=85^{\circ }$
(figure 4), in the region
$-0.5<\unicode[STIX]{x1D6FE}<0$
, and
$0<\unicode[STIX]{x1D714}_{r}<1$
. The colour plot is in logarithmic scale and denotes the value of
$\det (\unicode[STIX]{x1D63F})$
. One can visually identify the eigenvalues as the points where the concentric curves are centred. In order to measure quantitatively the error in the pseudospectra, we proceed as follows. We compute, for a fixed value of
$k=|\boldsymbol{k}|$
, the function
$f(\unicode[STIX]{x1D714})=\det (\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D714}))$
on a
$500\times 500$
grid covering the domain
$0\leqslant \unicode[STIX]{x1D714}_{r}\leqslant 5k$
and
$-0.5k\leqslant \unicode[STIX]{x1D6FE}\leqslant 0$
. Because
$f$
can vary by several orders of magnitude within this domain, it is not convenient to calculate an error directly defined on
$f$
. Moreover, we seek to define an error that emphasizes the differences only in regions where
$f(\unicode[STIX]{x1D714}_{VM})$
is small, disregarding regions where
$f$
is large. Therefore, we map
$f(\unicode[STIX]{x1D714})$
into a new function
$\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714})$
through a sigmoid:
$$\begin{eqnarray}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714})=1+\frac{9f}{1+f}.\end{eqnarray}$$
The new function
$\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714})$
is defined between 1 (for
$f=0$
) and 10 (for
$f\rightarrow \infty$
). Finally, we define the error
$\unicode[STIX]{x1D73A}_{D}$
as:
$$\begin{eqnarray}\unicode[STIX]{x1D73A}_{D}=\left\Vert \frac{(\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714})-\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714}_{VM}))}{\unicode[STIX]{x1D719}(\unicode[STIX]{x1D714}_{VM})}\right\Vert _{2},\end{eqnarray}$$
where the L-2 norm is normalized on the number of grid points.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=30^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=60^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Dispersion relation for ion-acoustic waves for angles of propagation
$\unicode[STIX]{x1D703}=80^{\circ },85^{\circ },89^{\circ }$
(curves from left to right). Blue, red and green lines are for fully kinetic, hybrid, and gyrokinetic models, respectively. (a) Real frequency
$\unicode[STIX]{x1D714}_{r}$
; (b) Damping rate
$\unicode[STIX]{x1D6FE}$
.
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis. All axes are in logarithmic scale.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=80^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=85^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
Errors for ion-acoustic wave with
$\unicode[STIX]{x1D703}=89^{\circ }$
, as functions of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Red and green curves are for hybrid and gyrokinetic models, respectively. Solid lines denote
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and dashed lines denote
$\unicode[STIX]{x1D73A}_{wave}$
.
2.2 Range of validity of the gyrokinetic model
As we mentioned, one important assumption in the derivation of the gyrokinetic model is that
$k_{||}\ll k_{\bot }$
. This requirement restricts dramatically the range of validity of the model. Following the philosophy outlined in the Introduction, we will compare the three models also in the non-quasi-perpendicular propagation regime, that is, in a range where the gyrokinetic model is not expected to be valid. Figure 5 shows the range of validity of the gyrokinetic model. The symbols denote the angle of propagation
$\unicode[STIX]{x1D703}$
for a given value of perpendicular wavevector
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
that results in a relative error on the kinetic Alfvén wave real frequency of
$10\,\%$
(squares) and
$30\,\%$
(diamonds). One can appreciate that, for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}>3$
, the propagation cone is restricted to angles very close to
$90^{\circ }$
. An alternative way of appreciating this result is shown in figure 6, which shows the dispersion relation for the Alfvén/ion-cyclotron wave. When plotted as
$\unicode[STIX]{x1D714}_{r}/k_{||}v_{A}$
as a function of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
, the dispersion relation becomes independent of angle, in the gyrokinetic model (shown in blue). The red, yellow and purple lines denote the fully kinetic solution for angles of propagation
$\unicode[STIX]{x1D703}=80^{\circ },85^{\circ },89^{\circ }$
. The gyrokinetic model has been suggested as a valid model for the study of solar wind turbulence all the way to electron scales (Howes et al.
Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011) on the basis that some observational data suggest that the propagation angle in the solar wind averages at approximately
$88^{\circ }$
and across the examined spectra never falls below
$60^{\circ }$
(Narita et al.
Reference Narita, Glassmeier, Sahraoui and Goldstein2010; Sahraoui et al.
Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010). It is not within the scope of this paper to discuss the issue of the observed wave spectrum anisotropy in the solar wind. However, it is fair to mention that the problem is still open, and some recent works have emphasized the importance of quasi-parallel wave propagation (Osman & Horbury Reference Osman and Horbury2009; Kiyani et al.
Reference Kiyani, Chapman, Sahraoui, Hnat, Fauvarque and Khotyaintsev2013; Lion, Alexandrova & Zaslavsky Reference Lion, Alexandrova and Zaslavsky2016). Furthermore, the development of wavevector anisotropy is an open research question, for which linear properties over a range of propagation angles will be important. It is in this spirit that, in the following, we present a comparison of the models across a wide range of propagation angles.
3 Results
Historically, different communities in plasma physics have adopted different nomenclatures for linear waves. Here we follow the nomenclature adopted by Gary (Reference Gary2005). We will refer to the Alfvén/ion-cyclotron mode as the electromagnetic mode that is left-handed polarized for parallel propagation. The ion-acoustic mode is the one that connects to the electrostatic solution for parallel propagation, and the fast (magnetosonic) waves is the mode that connects to the right-handed polarized mode for parallel propagation. In this section we first discuss the dispersion relations of each one of these three modes, and then we compare the models comparing the pseudospectra in the quasi-perpendicular regime.
3.1 Alfvén/ion-cyclotron (AIC) mode
Figure 7 shows the dispersion relation for the AIC mode, for angles of propagation
$\unicode[STIX]{x1D703}=10^{\circ },30^{\circ },60^{\circ }$
. Blue, red and green curves are for fully kinetic, hybrid and gyrokinetic models, respectively (the same notation will be used for all following figures). The real frequency
$\unicode[STIX]{x1D714}_{r}$
is shown in (a) and the damping rate in (b). In both panels
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis (logarithmic scale is used on all axes). The three different propagation angles can be identified from the starting values of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. Since all the curves are plotted from
$k\unicode[STIX]{x1D70C}_{i}=0.1$
, the initial perpendicular wavevector increases with increasing
$\unicode[STIX]{x1D703}$
. For
$\unicode[STIX]{x1D703}=10^{\circ }$
and
$30^{\circ }$
, the hybrid and fully kinetic solutions are almost indistinguishable. In particular, the saturation of wave frequency is well represented. The corresponding damping rates present a mismatch at low
$k$
and agree well at large
$k$
. A peculiar feature of the
$\unicode[STIX]{x1D703}=60^{\circ }$
solution is that a branch bifurcation happens around
$k_{\bot }\unicode[STIX]{x1D70C}_{i}=0.6$
for the hybrid solution. The gyrokinetic solution does not capture the correct frequency saturation, for any angle of propagation, as expected. Consequently, the frequencies of the gyrokinetic solutions become much higher than the correct values, even for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}<1$
. The gyrokinetic solutions are generally less damped than the correct fully kinetic solutions (b). The errors for the AIC mode solutions are shown in figures 8–10, for the cases
$\unicode[STIX]{x1D703}=10^{\circ }$
,
$30^{\circ }$
,
$60^{\circ }$
, respectively. Here, the hybrid model is shown in red and the gyrokinetic in green. Solid and dashed lines denote the errors
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and
$\unicode[STIX]{x1D73A}_{wave}$
, respectively. One can notice that the errors for the hybrid solutions are lower than the corresponding errors for the gyrokinetic model for
$\unicode[STIX]{x1D703}=10^{\circ }$
and
$30^{\circ }$
. However, which model performs better for
$\unicode[STIX]{x1D703}=60^{\circ }$
depends on
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
. The fact that the gyrokinetic errors in figure 10 for
$\unicode[STIX]{x1D703}=60^{\circ }$
do not grow very large, and almost saturate for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}>1$
despite the fact that the dispersion relation shows a large mismatch with the fully kinetic solution, is due to the fact that the mode becomes heavily damped, and therefore the contribution of the real frequency mismatch becomes negligible in the errors.
Moving to quasi-perpendicular angles, figure 11 shows the dispersion relation for AIC for
$\unicode[STIX]{x1D703}=80^{\circ },85^{\circ },89^{\circ }$
(from left to right). Interestingly, even at such a large propagation angle, the AIC mode resonates with ions and presents a frequency saturation. Once again the gyrokinetic solution is quite faithful (now both in frequency and damping rate) up to the wavenumbers where the dispersion relation saturates. In this regime, the hybrid solution presents larger errors, in particular close to the elbow of the curves where the damping rates depart from the Vlasov–Maxwell values (b). Moreover, in the
$\unicode[STIX]{x1D703}=89^{\circ }$
case, the saturated value of the hybrid solution is not correct (being larger than
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}_{i}=1$
). The respective errors are plotted in figures 12–14. It is clear how the gyrokinetic solution gets better and better, with respect to the hybrid, with increasing angle of propagation (only for wavevectors smaller than the value where dispersive effects take place).
3.2 Ion-acoustic (IA) mode
The dispersion relation for the IA mode is shown in figure 15, with the same format as in the previous section. One can distinguish the following features. Regarding the real frequency (a), there is again a good agreement between hybrid and fully kinetic; the gyrokinetic frequency also presents a similar agreement, but only for small
$k$
values. Damping rates present larger mismatches. Notably, the dispersion relation for the case
$\unicode[STIX]{x1D703}=60^{\circ }$
presents an ‘elbow’ for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$
that is captured by the hybrid model. We show the errors
$\unicode[STIX]{x1D73A}_{\unicode[STIX]{x1D714}}$
and
$\unicode[STIX]{x1D73A}_{wave}$
for
$\unicode[STIX]{x1D703}=10^{\circ }$
,
$\unicode[STIX]{x1D703}=30^{\circ }$
,
$\unicode[STIX]{x1D703}=60^{\circ }$
respectively in figures 16–18 (same format as previous section). The errors are now generally smaller for gyrokinetic mode at small
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
, and increasing with wavevector, and is almost constant in value for the hybrid solutions. The dispersion relation plots for quasi-perpendicular angles are shown in figure 19. Similarly to the AIC case, the gyrokinetic model reproduces more faithfully the Vlasov–Maxwell solutions, until dispersive effects take place. A small but noticeable mismatch in the sound speed can be noted for the hybrid solution, that is shifted vertically. The only curve where the gyrokinetic and the Vlasov–Maxwell solution agree for the whole range of wavevectors considered is for
$\unicode[STIX]{x1D703}=89^{\circ }$
. The corresponding errors are shown in figures 20–22. It is interesting to notice how the relative accuracy of the two models strongly depend on the wavevector. For instance, at
$\unicode[STIX]{x1D703}=80^{\circ }$
the gyrokinetic solution is approximately two order of magnitudes more accurate than the hybrid for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}=0.1$
, but starts to become less accurate for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}>2$
. Also, it is interesting that the hybrid solution errors decrease with increasing wavevector. That is, the model is more accurate in the dispersive regions of the curves.
Dispersion relation for fast waves for angles of propagation
$\unicode[STIX]{x1D703}=10^{\circ },30^{\circ },60^{\circ }$
(curves from left to right). Blue, red and green lines are for fully kinetic, hybrid and two-fluid models, respectively. (a) Real frequency
$\unicode[STIX]{x1D714}_{r}$
in logarithmic scale; (b) Damping rate
$\unicode[STIX]{x1D6FE}$
in logarithmic scale.
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
is on the horizontal axis in logarithmic scale.
3.3 Fast mode
It is well known that the fast mode solution is factored out from the gyrokinetic model. In this sense, a comparison of the gyrokinetic with the other the two models is meaningless. In figure 23 we compare the dispersion relations of the fully kinetic, the hybrid and the two-fluid models for
$\unicode[STIX]{x1D703}=10^{\circ }$
,
$\unicode[STIX]{x1D703}=30^{\circ }$
,
$\unicode[STIX]{x1D703}=60^{\circ }$
. The dispersion relation for the latter follows the derivation of Zhao et al. (Reference Zhao, Voitenko, Yu, Lu and Wu2014). The three solutions are in a very good agreement for the real frequency (a), where a slight deviation with respect to the correct solution occurs only at large
$k_{\bot }$
. Interestingly, the two-fluid solution is a very good proxy for both hybrid and fully kinetic solutions, in this regime. (b) Shows the corresponding damping rate (the two-fluid model is not shown since its solutions are always real). Once again, as expected, the hybrid model underestimates the damping rates, although there is a certain regime, approximately for
$k\unicode[STIX]{x1D70C}_{i}\leqslant 1$
, where the red and blue curves have a good agreement.
Error
$\unicode[STIX]{x1D73A}_{D}$
as a function of
$k_{\bot }\unicode[STIX]{x1D70C}_{i}$
for angles of propagation
$\unicode[STIX]{x1D703}=60^{\circ },70^{\circ },80^{\circ },85^{\circ }$
(blue, red, yellow and magenta). Dashed and solid lines are for gyrokinetic and hybrid models, respectively.
Real (black) and imaginary (red) part of
$f(\unicode[STIX]{x1D714})=\det (\unicode[STIX]{x1D63F})$
for
$\unicode[STIX]{x1D703}=85^{\circ }$
,
$k\unicode[STIX]{x1D70C}_{i}=1$
. (a) Fully kinetic; (b) Hybrid; (c) Gyrokinetic model. The intersection between black and red curves identify normal modes.
3.4 Quasi-perpendicular propagation
In this section we perform the analysis of the error
$\unicode[STIX]{x1D73A}_{D}$
of (2.6). As we discussed above, in this regime the goal is not to compare the solutions for normal modes, but rather to provide a single measure of the error that captures the overall distance between the models in term of their pseudospectrum. This approach is justified by the non-normality of the linear operators, and by the fact that a multitude of modes appear in this regime, and it is not obvious why any of those should play a particular role in the plasma dynamics. Figure 24 shows the error
$\unicode[STIX]{x1D73A}_{D}$
for angles of propagation
$\unicode[STIX]{x1D703}=60^{\circ }$
,
$\unicode[STIX]{x1D703}=70^{\circ }$
,
$\unicode[STIX]{x1D703}=80^{\circ }$
and
$\unicode[STIX]{x1D703}=85^{\circ }$
(respectively in blue, red, yellow and magenta). The solid lines are for the hybrid model, and the dashed lines for the gyrokinetic. A very interesting trend can be noticed. The errors for the gyrokinetic model have a sharp increase, of almost an order of magnitude, around
$k_{\bot }\unicode[STIX]{x1D70C}_{i}=1$
, while the hybrid model error tends to increase for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}<1$
and to decrease for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}>1$
. Overall, the hybrid model performs much better than the gyrokinetic. This is not surprising because, as was already evident by inspection of figures 3 and 4, the gyrokinetic model simply misses many roots of the fully kinetic dispersion relation. It is important to realize that the modes that are not captured by this model are not ‘heavily damped’ with respect to the modes that are captured (but they have higher frequencies). One could argue that the error
$\unicode[STIX]{x1D73A}_{D}$
is defined and constructed in a rather arbitrary way. Although this is certainly true, we think that it is still very informative regarding how much a reduced model is able to reproduce the correct linear evolution of the system. As we mentioned, the ability of the gyrokinetic model to describe solar wind turbulence in the regime
$k_{\bot }\gg k_{||}$
has been stressed in many recent papers (Schekochihin et al.
Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; TenBarge, Howes & Dorland Reference TenBarge, Howes and Dorland2013). Although this might certainly be true, we believe that it is informative to know what is left out of the model in terms of normal modes. In figure 25 we compare the three models (fully kinetic to the left, hybrid in the middle, gyrokinetic to the right) for the case
$k\unicode[STIX]{x1D70C}_{i}=1$
and
$\unicode[STIX]{x1D703}=85^{\circ }$
. Here, we do not show the normal modes but the curves
$\text{Re}(f)=0$
(black) and
$\text{Im}(f)=0$
(red) in the domain
$-0.2<\unicode[STIX]{x1D6FE}<0$
and
$0<\unicode[STIX]{x1D714}_{r}<2.5$
. The normal modes can be inferred as the points where black and red curves intersect. This figure reinforces what is well known about the models: the gyrokinetic model factors out high-frequency modes (irrespective of their damping rate) and the hybrid model does not properly damp some of the modes (the ones that undergo electron Landau damping).
4 Conclusions
Reduced kinetic models are extensively used in plasma physics in lieu of the more complete Vlasov–Maxwell system. In this paper, we have compared the hybrid (ion kinetic and electron fluid), the gyrokinetic and the fully kinetic models, for a wide range of wavevectors and propagation angles, in linear theory. We have chosen one test case with parameters relevant to the study of solar wind turbulence. We have compared the dispersion relations of single normal modes, for angles of propagation
$\unicode[STIX]{x1D703}=10^{\circ },30^{\circ },60^{\circ },80^{\circ },85^{\circ },89^{\circ }$
and pseudospectra for
$\unicode[STIX]{x1D703}=60^{\circ },70^{\circ },80^{\circ },85^{\circ }$
. The latter method is justified by the multitude of modes that appear in the quasi-perpendicular case, which makes it non-trivial to identify which one plays a predominant role. It is important to remember that the gyrokinetic model has been employed also outside its nominal range of validity (
$k_{||}\ll k_{\bot }$
). We note that the Vlasov–Maxwell linear theory has recently been investigated by Sahraoui et al. (Reference Sahraoui, Belmont and Goldstein2012), for very oblique propagation angles. Our findings can be summarized as follows. The hybrid solutions present a reasonable agreement with the Vlasov–Maxwell system across all propagation angles, the largest disagreement being for
$\unicode[STIX]{x1D703}=89^{\circ }$
. Interestingly, the errors decrease in the region where the waves resonate with ions and the dispersion curves saturate. One major mismatch of the hybrid model happens for
$\unicode[STIX]{x1D703}=60^{\circ }$
, where a clear matching of the hybrid AIC and IA modes with the corresponding Vlasov–Maxwell modes is not possible due to branch switching. As expected, the gyrokinetic model presents larger errors for angles of propagation less than
$80^{\circ }$
, particularly for
$k_{\bot }\unicode[STIX]{x1D70C}_{i}>1$
. However, also for quasi-perpendicular propagation, the errors increase abruptly as soon as the waves become resonant with the ions. In conclusion, it appears that the hybrid model should be preferred for its ability to represent most of the normal modes, at least in the ion regime
$k_{\bot }\unicode[STIX]{x1D70C}_{i}<1$
and, to a certain extent, also in the electron regime. However, we notice that the spurious low value of damping rates of some hybrid solutions might result in numerical instabilities, due to the accumulation of energy that is not properly damped. In this respect, a promising numerical approach consists in devising models that factor out Langmuir oscillations while still retaining the fully kinetic nature of the solution (see, e.g. Degond, Deluzet & Savelief Reference Degond, Deluzet and Savelief2012; Chen & Chacon Reference Chen and Chacon2014; Tronci & Camporeale Reference Tronci and Camporeale2015).
