text:
\section{Introduction}

One of the long-pursued but yet-to-be discovered particles are Majorana fermions.
They are hypothetical elementary particles \cite{Elliott2015}, that satisfy the condition that their antiparticles are identical to their particles.
Although such particles have not been found in particle physics, recent theoretical studies surprisingly suggest that Majorana fermions can be generated as emergent excitations in condensed matter systems \cite{Read2000,Qi2011,Sato2016,Sato2017,Alicea2012,Beenakker2013,Aguado2017}.
Thus, the naturally arising questions are whether or not such exotic particles really exist in nature and what novel properties they have if they exist.

Fermi superfluids and superconductors are natural candidates for systems that host Majorana fermions.
In these systems, quasiparticles (QPs) are described as superpositions of a particle and a hole state, and therefore the QPs could be identical to their antiparticles if certain conditions are satisfied \cite{Read2000,Qi2011,Sato2016,Sato2017,Alicea2012,Beenakker2013,Aguado2017}.
In fact, Majorana fermions are predicted to appear as emergent excitations at a surface or in a vortex of topological superfluids/superconductors.
Majorana fermions have attracted considerable attention for the last decade not only because of fundamental interest in their properties but also for their potential application to fault-tolerant quantum computing, and much experimental effort has been devoted to the search for Majorana fermions  \cite{Sato2016,Sato2017,Beenakker2013}.

The principal systems that host Majorana excitations are odd-parity topological superfluids/superconductors.
The $p$-wave superfluid $^3$He is a prime candidate.
For superconductors, there are several candidates which are now conceived to have a nontrivial topology in the bulk order parameters \cite{Yonezawa2016}, such as Sr$_2$RuO$_4$ \cite{Maeno2012}, the B phase of UPt$_3$ \cite{Izawa2014,Goswami2015}, Cu$_x$Bi$_2$Se$_3$ \cite{Sasaki2015}, and Sn$_{1-x}$In$_x$Te \cite{Sasaki2015}.
In some of these superconductors, the presence of edge or surface Andreev bound states has been reported by tunneling spectroscopy measurements (Sr$_2$RuO$_4$ \cite{Kashiwaya2011}, Cu$_x$Bi$_2$Se$_3$ \cite{Sasaki2011}, and Sn$_{1-x}$In$_x$Te \cite{Sasaki2012}), suggesting topologically nontrivial states in the bulk.
However, the bulk order parameters have not yet been firmly established for these superconductors.
Other intensively studied systems is a nanowire-superconductor hybrid system, where a topological superconducting state is induced in the nanowire \cite{Oreg2010,Lutchyn2010}.
Zero-energy states have been observed at the ends of nanowires by tunneling spectroscopy \cite{Mourik2012,Das2012,Deng2012,Deng2016,Albrecht2016}, but there is an unsettled question as to whether they are trivial or nontrivial zero modes at present \cite{Lee2014}.

In order to make clear understanding of properties of Majorana fermions, experiments in well-elucidated topological systems are essential.
The superfluid $^3$He \cite{Vollhardt90,Dobbs00,Volovik2003} provides an ideal platform for exploring the nature of Majorana fermions because its bulk properties are well elucidated and clear experimental results are obtainable owing to the cleanliness of liquid $^3$He.
One of the superfluid phases, called the B phase ($^3$He-B), is recognized as a topological superfluid protected by time-reversal symmetry \cite{Mizushima2016,Schnyder2008,Chung2009}.
On the surface of $^3$He-B, Majorana fermions are expected to emerge as surface Andreev bound states \cite{Schnyder2008,Chung2009,Qi2009,Volovik2009_1}.

The surface Andreev bound states in $^3$He-B have been observed in the last decade by transverse acoustic impedance and other methods \cite{Aoki2005,Okuda2012,Murakawa2009,Murakawa2011,Choi2006}.
Their observations are very important because they verify that the nontrivial topology of $^3$He-B generates such bound states.
In these studies, however, the employed surfaces have more or less microscopic roughness,  which prevents understanding of detailed properties of surface Andreev bound states.
At a smooth surface called a specular surface, on the other hand, surface Andreev bound states with a well-defined energy spectrum are formed, and exotic features of Majorana fermions should be unambiguously elucidated.
A specular surface is realized at a free surface (see Sect. \ref{sect:specularity}).
In order to detect Majorana fermions formed at a free surface, Ikegami {\it et al.} measured the mobility of ions trapped below a free surface of $^3$He-B in 2013 \cite{Ikegami2013-2}.
Recently, the experimental mobility was quantitatively reproduced by a calculation by Tsutsumi that included contributions from surface Andreev bound states \cite{Tsutsumi2017}.
This quantitative agreement unambiguously demonstrates the experimental detection of surface Andreev bound states at a free surface.
In this review, we describe such recent progress with the aim of providing a clear view of the detection of surface Andreev bound states from the mobility of ions.

\section{Majorana Fermions Formed at a Surface of $^3$He-B}
\label{Majorana_3He}

The superfluid $^3$He is a $p$-wave pairing state with orbital angular momentum 1 and spin angular momentum 1.
In $^3$He-B, the condensate is composed of an admixture of Cooper pairs with three substates of the orbital angular momentum ($L_z=$ $-$1, 0, 1) and three substates of the spin ($S_z=$ $-$1, 0, 1), and
the orbital angular momentum $\mbox{\boldmath $L$}$ and the spin angular momentum $\mbox{\boldmath $S$}$ couple with each other such that the twisted total angular momentum becomes zero [$\mbox{\boldmath $L$}+R(\mbox{\boldmath $\hat n$},\theta)^{-1}\mbox{\boldmath $S$}=0$, where $R(\mbox{\boldmath $\hat n$},\theta)$ is the rotation matrix].
The order parameter $\mbox{\boldmath $d$}$ is then described as \cite{Vollhardt90,Dobbs00}
\begin{equation}
\mbox{\boldmath $d$} (\mbox{\boldmath $\hat k$}) = \Delta_B R(\mbox{\boldmath $\hat n$},\theta) \mbox{\boldmath $\hat k$} ,
\label{eq:OP}
\end{equation}
where $\Delta_B$ is the energy gap of $^3$He-B and $\mbox{\boldmath $\hat k$}$ is a unit vector in momentum space.
This state breaks the relative spin-orbit rotation symmetry because of the coupling of $\mbox{\boldmath $L$}$ and $\mbox{\boldmath $S$}$.
This state is however invariant under time inversion operation $\hat T$ because $\hat T \mbox{\boldmath $d$} (\mbox{\boldmath $\hat k$}) = -\mbox{\boldmath $d$}^* (-\mbox{\boldmath $\hat k$})= \mbox{\boldmath $d$} (\mbox{\boldmath $\hat k$})$.

One of the notable phenomena expected in $^3$He-B is the formation of surface Andreev bound states in the region within the coherence length from the surface.
In 1981, Buchholtz and Zwicknagl first pointed out that $s$-wave and $p$-wave superconductors respond to a surface differently and showed that Andreev bound states are formed at the surface of $^3$He-B \cite{Buchholtz1981}.
Since then, a number of theoretical investigations not only for a specular surface but also for surfaces in the presence of microscopic roughnesses have been carried out for $p$-wave \cite{Nagato1998,Zhang1987,Zhang1988,Nagai2008} and $d$-wave pairing states \cite{Hu1994,Matsumoto1995,Kashiwaya1995}.
Basically, surface Andreev bound states are formed as a consequence of the interference of an incident wavefunction and a reflected wavefunction of a QP \cite{Ohashi1996}.

The surface Andreev bound states of  $^3$He-B are intimately related to the non-trivial topology of the bulk $^3$He-B order parameter.
In the superfluid $^3$He, investigations of the topological aspects  have already begun as early as the 1980's \cite{Salomaa1988}.
After that, a number of unusual behaviors have been uncovered by using the concept of topology \cite{Volovik2003}, especially in the last decade, along with the drastic progress in the topological physics in condensed matter systems \cite{Qi2011,Sato2016,Sato2017,Beenakker2013,Mizushima2016,Hasan2010,Ando2013}.
Now, $^3$He-B is recognized as a topological state protected by time-reversal symmetry \cite{Mizushima2016,Schnyder2008,Chung2009,Mizushima2015}.
In general, the nontrivial topology of a bulk gives rise to formation of gapless bound states on its surface owing to the so-called bulk-surface correspondence.
In $^3$He-B, gapless surface Andreev bound states are formed as a result of the broken symmetry of the bulk $^3$He-B order parameter.
Because of equal-weight superpositions of a particle and a hole state, surface Andreev bound states satisfy the Majorana condition as shown below \cite{Schnyder2008,Chung2009}.

To see basic properties of surface Andreev bound states, we consider solutions of the Bogoliubov-de Gennes equation near a surface.
We take the $z$-axis normal to the surface of $^3$He-B which fills the region at $z>0$.
We approximate the perpendicular and the parallel components of the order parameter as $\Delta_\bot(z)=\Delta_B \tanh (\kappa z)$ and $\Delta_\parallel (z)=\Delta_B$ near the surface respectively ($\kappa  = \Delta _B /\hbar v_F$) \cite{Tsutsumi2012-2}, which provides a good approximation for the suppression of the order parameter near the surface (see Appendix in Ref. \cite{Tsutsumi2012-2}).
The field operator for the surface Andreev bound states is then given by \cite{Tsutsumi2017}
\begin{equation}
\begin{array}{rr}
 \left[ {\begin{array}{*{20}c}
   {\psi _ \to  ({\bf{r}})}  \\
   {\psi _ \leftarrow  ({\bf{r}})}  \\
   {\psi _ \to ^\dag  ({\bf{r}})}  \\
   {\psi _ \leftarrow ^\dag  ({\bf{r}})}  \\
\end{array}} \right] & = \sum \limits_{\bf{k}} {A_k \frac{1}{\cosh (\kappa z)} \sin (k_ \bot  z)\left( {\hat \gamma _{{\bf{k}} } e^{i{\bf{k}}_\parallel \cdot {\bf{r}} }  + \hat \gamma _{{\bf{k}} }^\dag  e^{ - i{\bf{k}}_\parallel  \cdot {\bf{r}} } } \right)}  \\
  & \times \left[ {\begin{array}{*{20}c}
   {\cos \left[ {(\phi _k  + \pi /2)/2} \right]}  \\
   {\sin \left[ {(\phi _k  + \pi /2)/2} \right]}  \\
   {\cos \left[ {(\phi _k  + \pi /2)/2} \right]}  \\
   {\sin \left[ {(\phi _k  + \pi /2)/2} \right]}  \\
\end{array}} \right] . \\
 \end{array}
\label{eq_SABS}
\end{equation}
Here $\hat \gamma _{{\bf{k}} }^\dag$ ($\hat \gamma _{\bf{k}}$) is a Fermion creation (annihilation) operator, ${\bf{k}}=({\bf{k}}_\parallel,k_\bot)$, ${\bf{k}}_\parallel$ is the momentum parallel to the surface,  $k_\bot^2=k_F^2-{\bf{k}}_\parallel ^2$,  $A_k$ is a normalization constant, and $\to $ ($\leftarrow$) represents the spin $\sigma$ with a quantization axis taken in the $x$-direction ($\hbar$ is Planck's constant divided by 2$\pi$, $v_F$ is the Fermi velocity, and $k_F$ is the Fermi wavenumber).
These states localize in the vicinity of the surface with the length scale of $\kappa^{-1}=3.6\xi_0$, where $\xi_0=\hbar v_F/2\pi k_BT_c \sim$ 80 nm is the coherence length ($k_B$ is Boltzmann's constant and $T_c$ is the superfluid transition temperature).
The modes described in Eq. (\ref{eq_SABS}) have a gapless linear dispersion relation with energy $E$ given by $E= \hbar c \left| \bf{k}_\parallel \right|$ with $c=\Delta_B/(\hbar k_F)$.
This dispersion relation is called Majorana cone.
As shown in Refs. \cite{Chung2009,Mizushima2015,Mizushima2016}, the surface Andreev bound states are described by the massless two-dimensional (2D) Dirac Hamiltonian.
From the linear dispersion relation, the local density of states of the surface Andreev states is obtained as $N_s(E,z)= \frac{\pi}{4} N_0 \frac{E}{\Delta_B} \frac{1}{[\cosh (\kappa z)]^{2}}$, where $N_0$ is the density of states at Fermi energy in normal $^3$He \cite{Tsutsumi2012-2}.

The remarkable feature of surface Andreev bound states is that they are described as Majorana fermions.
In fact, they behave as relativistic QPs described by the massless 2D Dirac Hamiltonian \cite{Chung2009,Mizushima2015,Mizushima2016} and satisfy the Majorana condition, $\psi _ \sigma ({\bf{r}}) = \psi _  \sigma ^\dag  ({\bf{r}})$, as seen in Eq. (\ref{eq_SABS}).

\begin{figure}
\begin{center}
\includegraphics[width=0.7\linewidth,keepaspectratio]{surface_condition.eps}
\end{center}
\caption{\label{surface_condition} (Color online) Schematic pictures of (a) specular reflection and (b) diffusive reflection for a QP in normal state.
}
\end{figure}

The properties of surface Andreev bound states are strongly affected by nature of the surface.
In the case of an ideal surface called specular surface [Fig.  \ref{surface_condition}(a)], a QP in normal state is reflected in the angle same to the incident angle.
In fact, the wavefunction described by Eq. (\ref{eq_SABS}) is the solution for the specular surface;
 the modes in Eq. (\ref{eq_SABS}) are eigenfunction of the momentum parallel to the surface ${\bf{k}}_{\parallel}$.
(Note that the momentum parallel the surface is conserved in the specular reflection.)
Such a specular surface is realized at a free surface of liquid $^3$He, as shown in Sect. \ref{sect:specularity}.
Another important surface often used in experimental studies is a wall of a solid with which liquid $^3$He is in contact [Fig.  \ref{surface_condition}(b)].
In this case, irregularities in atomic scale ($\sim k_F^{-1}$) existing on the surface make QPs reflected in random directions.
In this diffusive reflection, ${\bf{k}}_{\parallel}$ is not conserved, which makes the spectra of surface Andreev bound states to be broadened and shifted significantly \cite{Nagato1998,Nagai2008}.
Therefore, surface Andreev bound states in this case no more have a well-defined Majorana cone.
Note that, in the diffusive reflection, the momentum  on the order of $p_F (=\hbar k_F)$ parallel to the surface is exchanged between the wall and the QP.

Experimentally, surface Andreev bound states formed on solid surfaces have been observed so far.
After the first detection by transverse acoustic impedance measurement in 2005 \cite{Aoki2005}, they have been observed not only by this method \cite{Okuda2012,Murakawa2009,Murakawa2011} but also by specific heat \cite{Choi2006}, by the attenuation of transverse sound \cite{Davis2008}, and potentially by damping and critical velocity of a micromechanical oscillator \cite{Zheng2016,Zheng2017}.
In particular, a series of studies using transverse acoustic impedance \cite{Okuda2012,Murakawa2009,Murakawa2011} showed, although indirectly, that the dispersion relation of surface Andreev bound states approaches to the Majorana cone by changing the surface quality from diffusive to partially specular by coating the surface with $^4$He layers;
The experiments themselves cold not detect the dispersion relation, but the observed acoustic impedance was consistently explained by the theory in which the dispersion relation is assumed to approach to the Majorana cone as increasing specularity.
However, even if the specularity is enhanced by the $^4$He coating, the surface cannot be perfectly specular for a solid surface, which prevents us from quantitative understanding of nature of surface Andreev bound states.
On the other hand, at a specular surface, where surface Andreev bound states should be Majorana fermions having a well-defined linear dispersion, one can reach unambiguous understanding of their properties.

\section{Specular Nature of Free Surface of Liquid $^3$He}
\label{sect:specularity}

As mentioned in Sect. \ref{Majorana_3He}, Majorana surface Andreev bound states have a linear dispersion relation only at a specular surface.
A specular surface can be realized at a free surface of liquid $^3$He.
The free surface is an interface in contact with $^3$He vapor (essentially a vacuum in the sub-mK region), and there is no impurity or disturbance at the surface, except for thermally excited surface waves.
The thermal surface waves, however, do not degrade the specular nature of reflection.
To make the reflection diffusive, irregularity of the surface at a length scale on the order of the Fermi wavelength $\lambda _F =2\pi/ k_F$ ($=$ 0.82 nm)  is necessary, but the wavelength of thermally excited surface waves is, at 1 mK for example, $\lambda _{\rm th} =2\pi (\alpha_3/\rho_3)^{1/3}(\hbar /k_BT)^{2/3} \sim$ 0.3 $\mu$m, much longer than $\lambda _F$.
(Here we used the dispersion relation of capillary waves $\omega=(\alpha_3/\rho_3)^{1/2}k^{3/2}$ with frequency $\omega$ and wavenumber $k$, where $\alpha_3$ and $\rho_3$ are the surface tension and density of liquid $^3$He, respectively.)
Therefore, the free surface can be regarded as smooth at the atomic scale even in the presence of thermally excited surface waves, suggesting that QPs are reflected specularly.

The specular nature of a free surface of liquid $^3$He has actually been demonstrated via the mobility $\mu$ of a Wigner crystal of electrons floating above a free surface \cite{Book2DE,BookMonarkha}.
In the Wigner crystal phase, which is formed at sufficiently low temperatures (the transition temperature is 0.32 K for an electron density of 1$\times$10$^{12}$ m$^{-2}$ for example), the electrons exert a spatially modulated pressure on the surface, generating a deformation of the surface commensurate with the Wigner crystal called a dimple lattice.
Although the depth of the dimple  $\delta$ is only $\sim$ 0.1 \AA , the dimple lattice significantly affects the mobility of elections parallel to the surface.
In particular, the motion of the Wigner crystal involves the motion of the dimple lattice, which is impeded by a drag force acting on the dimple lattice.
In the ballistic regime realized at low temperatures ($T \lesssim$ 20 mK), the drag force is caused by the reflection of ballistic QPs at the surface.
In this ballistic regime, experimental studies have revealed that $\mu$ in normal $^3$He is independent of $T$ \cite{Shirahama1997,Monarkha1998}, and $\mu$ in $^3$He-B and $^3$He-A rapidly increases with decreasing $T$ \cite{Shirahama1997,Monarkha1998,Ikegami2006}.
These experimental mobilities are quantitatively explained by the theoretical model of specular reflection \cite{Monarkha1997}, demonstrating the specular reflection of QPs at the free surface.

The quantitative agreement of the mobility with the specular reflection model indicates that the fraction of diffusively reflected QPs should be very small (less than one in 10$^{10}$ of the incident QPs).
This can be understood by noting that even such a small fraction of diffusively-reflected QPs can greatly reduce the mobility.
To see this, we consider a drag force $dF_D$ acting on a surface element $dS$ of a dimple lattice moving horizontally with velocity $V_0$.
The drag force in the case of specular reflection is given as
\begin{equation}
 dF_D^{(s)}  \sim n_3 p_F V_0 (\delta /a)^2 dS ,
\label{eq:specular}
\end{equation}
while that for diffusive reflection is
\begin{equation}
dF_D^{(d)} \sim n_3 p_F V_0 dS ,
\label{eq:diffuive}
\end{equation}
where $a$ is the period of the dimple lattice ($\sim$ 1 $\mu$m) and $n_3$ is the number density of liquid $^3$He \cite{Monarkha2006} (see also the argument described in Ref. \cite{Ikegami2017}).
Equations (\ref{eq:specular}) and (\ref{eq:diffuive}) suggest that the diffusive reflection generates a drag force that is larger by a factor of $(a/ \delta)^2$ ($\sim$ 10$^{10}$).
Therefore, the mobility should be significantly suppressed even if a tiny fraction of QPs on the order of 10$^{-10}$ undergo conventional diffusive reflection.

\section{Ions and Experimental Details}
\subsection{Ions in Liquid $^3$He}
\label{sect:ion}

Ions in liquid $^3$He are charged objects at the mesoscopic scale, which have been used to investigate the microscopic properties of liquid $^3$He.
There are two species of ions: a negatively charged ion called an electron bubble and a positively charged ion called a snowball \cite{Fetter1976}.
In this article, we only consider the negative ion.
The negative ion is the state of an electron self-trapped in a small spherical void with a radius of about 1 nm \cite{Fetter1976}.
The negative ion can be generated by injecting an electron into liquid helium, and the void is created as a result of the Pauli exclusion principle between the electron and electrons in the shells of $^3$He atoms.
The size of the void is determined by the balance of the zero-point energy of the confined electron and the surface tension of the void at zero pressure.
The void can be accurately regarded as a hard-sphere potential for QPs \cite{Shevtsov2017} because QPs cannot penetrate inside the void due to a huge repulsive barrier.
The hard-core radius $R$ is estimated to be $R=11.17k_F^{-1}$ ($k_F^{-1}=$ 0.13 nm) \cite{Shevtsov2016} from the low-temperature mobility in normal $^3$He at zero pressure obtained by Ikegami {\it et al.} \cite{Ikegami2013-2,Ikegami2015}.

Transport measurements of the negative ions have been used to uncover unusual aspects of the superfluid $^3$He.
In bulk $^3$He-B, the mobility has revealed that the $p$-wave coherence of QPs causes an unusual reduction in the transport cross section of the ions \cite{Baym1977,Baym1979,Ahonen1978}.
In $^3$He-A, transport of the ions was found to exhibit an anomalous Hall effect \cite{Ikegami2015,Ikegami2013-1} as a result of the skew scattering of QPs by the ions \cite{Shevtsov2016}, which directly demonstrated the breaking of time-reversal symmetry.
(We note that Refs.\cite{Salmelin1989,Salmelin1990} pointed out that the chiral order parameter of $^3$He-A could allow for transverse Hall term in mobility tensor.
However, the formulation described there gave zero Hall term \cite{Shevtsov2016}.)
In this article, we show that the transport properties of the negative ion are also useful for detecting surface Andreev bound states.

\subsection{Ions Trapped below Free Surface}
\label{sect:ion_below_surface}

Ions can be trapped below the free surface in a potential well $U(z) $ produced by the combination of the image charge of an ion and an external vertical electric field $E_{ \bot }$ \cite{Poitrenaud1972} [Fig. \ref{trapped_ion}(b)]:
\begin{equation}
U(z)  = \frac{{1}}{{16\pi \varepsilon \varepsilon _0 }}\left( {\frac{{\varepsilon  - 1}}{{\varepsilon  + 1}}} \right)\frac{e^2 }{z} + eE_{\bot}  z ,
\label{potential}
\end{equation}
where $z$ is the distance from the surface, $\varepsilon _0$ is the vacuum permittivity, and $\varepsilon$ is the relative permittivity of liquid $^3$He.
In the potential well, the ions are located at the minimum of $U(z)$ at $z=d$.
The position of the minimum can be adjusted in the range of 20 $< d <$ 60 nm by tuning $E_{ \bot }$ under typical experimental conditions.
This tunability of the depth provides an interesting opportunity for detecting surface Andreev bound states:
The amplitude of the wavefunction of surface Andreev bound states varies over $\xi_0$ ($\sim$ 80 nm) [see  Fig. \ref{trapped_ion}(a)], which should cause a change in the collision frequency of the ions with surface Andreev bound states if the depth is changed.

\begin{figure}
\begin{center}
\includegraphics[width=0.6\linewidth,keepaspectratio]{trapped_ion.eps}
\end{center}
\caption{\label{trapped_ion} (Color online) (a) Surface Andreev bound states and ions trapped below the free surface. The surface Andreev bound states are formed in the region within $\sim \xi_0 $ ($\sim$ 80 nm) from the surface, and the ions are trapped at a depth that is tunable in the range of 20--60 nm. (b) Trap potential for an ion as a function of distance $z$ from the surface. The curves represent the potential [Eq. (\ref{potential})] at $E_{ \bot }$=2, 5, and 10 kV/m. The arrows indicate the minimum in the potential, at which the ions are trapped.
}
\end{figure}

\subsection{Experimental Details of Mobility Measurements of Ions below Free Surface}
\label{sect:Exp_ion_surface}

In the experiments described in Ref. \cite{Ikegami2013-2}, the mobility $\mu$ of ions trapped below a liquid $^3$He surface was measured with a Corbino disk composed of two concentric electrodes attached to the ceiling of an experimental cell.
A bottom circular electrode, which was located 3.0 mm below the Corbino disk, was used to provide a vertical electric field $E_{ \bot }$ to control the depth of the trapped ions.
The free surface was set around the midpoint between the Corbino disk and the bottom electrode.
Negative ions were produced by field emission from a tungsten tip immersed in the liquid $^3$He, and were trapped below the surface.
The trapped ions formed a 2D ion sheet with a density of 1.5$\times$10$^{11}$ m$^{-2}$.
The mobility was measured using the Sommer--Tanner technique.
An ac voltage with a frequency of 0.1--12 Hz was applied on one of the Corbino electrodes, and the induced current was recorded by the other electrode.
The mobility was deduced from the longitudinal conductivity obtained from the output current.
The velocity of the ions was in the range of 10$^{-4}-$10$^{-3}$ m/s, sufficiently lower than the pair-breaking velocity ($\sim$10$^{-2}$ m/s).

The liquid $^3$He was cooled to $\sim$ 250 $\mu$K by a heat exchanger made of packed silver and platinum powders.
The temperature $T$ was measured by a platinum NMR thermometer mounted on a nuclear stage above 500 $\mu$K.
Below 500 $\mu$K, the temperature of liquid $^3$He was directly determined by the density of thermally excited QPs in bulk $^3$He-B measured with a vibrating wire immersed in the liquid.
To operate the vibrating wire, a magnetic field of 30 mT was applied perpendicular to the surface.
This field should open a Zeeman gap of $\sim$ 20 $\mu$K for surface Andreev bound states, but it does not affect the mobility because the temperature range for the measurements is much higher.
We also note that the data taken at 0 mT showed the same temperature dependence as described in Ref. \cite{Ikegami2013-2}.
For more details of the experiments, see Ref. \cite{Ikegami2013-2}.

\section{Detection of Surface Andreev Bound States}

\begin{figure}
\begin{center}
\includegraphics[width=0.7\linewidth,keepaspectratio]{3He-B_mobility.eps}
\end{center}
\caption{\label{3He-B_mobility} (Color online)
Mobility of negative ions trapped below the surface of $^3$He-B.
Here the mobility is normalized by the value $\mu_c=$1.8$\times$10$^{-6}$ m$^2$/Vs at the transition temperature  $T_c$(= 0.93 mK).
Inset: $\mu/ \mu_c$ as a function of $T_c/T$.
The experimental data were taken at four different depths, and the mobilities measured at these depths show the same temperature dependence.
The blue curve represents the theoretical mobility in bulk $^3$He-B calculated by Baym {\it et al} \cite{Baym1977,Baym1979}.
The red, orange, and light blue curves are, respectively, the theoretical mobility calculated at depths of 58, 40, and 21 nm by including the contribution from surface Andreev bound states \cite{Tsutsumi2017}.
Because of the lack of depth dependence, these theoretical curves fall on a single curve.
The theoretical curves are for $k_FR=$ 11.17.
The experimental data are taken from Ref. \cite{Ikegami2013-2} and the theoretical curves are reproduced from Ref. \cite{Tsutsumi2017}.
The shaded area corresponds to the contribution from surface Andreev bound states.
}
\end{figure}

The basic idea of the experiment to detect surface Andreev bound states using the mobility of the ions is as follows \cite{Ikegami2013-2}.
If the mobility is measured at a depth $d< \xi_0$, the mobility could be affected by scatterings with surface Andreev bound states, which are formed within a distance of $\xi_0$ from the surface [Fig. \ref{trapped_ion}(a)].
At sufficiently low temperatures, the mobility could be dominated by scattering of thermally excited surface Andreev bound states because the density of thermally excited QPs in bulk rapidly decreases with decreasing temperature.
If so, the mobility should be more suppressed when it is measured closer to the surface because of the increased scattering of surface Andreev bound states owing to the larger amplitude of the wavefunction.

The mobility of negative ions measured below the free surface of $^3$He-B is shown in Fig. \ref{3He-B_mobility} \cite{Ikegami2013-2}.
The mobility considerably increases with decreasing temperature.
In the figure, the mobilities measured at four different depths in the range of 21--58 nm are exhibited.
Interestingly, there is no depth dependence of the mobility within $\pm$1\% in this depth range.

The mobility should be affected not only by thermally excited surface Andreev bound states but also by thermally excited bulk QPs.
The mobility limited by bulk QPs was calculated by Baym {\it et al.} \cite{Baym1977,Baym1979}, which is shown by the blue line in Fig. \ref{3He-B_mobility}.
It increases rapidly with decreasing temperature because of the exponential reduction of the bulk QP density $n_{\rm QP}$($\sim e^{-\Delta_B /k_BT}$) and the reduction of the transport cross section of the ions for  QPs with low energies \cite{Baym1977,Baym1979}.
The mobility measured below the surface appears to be suppressed from the theoretical mobility in the bulk at $T/T_c <$ 0.5.
(We note that there has been no measurement of mobility in the bulk at $T/T_c \lesssim$ 0.8 at low pressures \cite{Ahonen1978,Ahonen1976}, thus we compare the data with the theoretical mobility in the bulk.)
This suppression implies that the mobility is affected by surface Andreev surface states at low temperatures.

If there are some contributions from surface Andreev bound states, the mobility is expected to be more suppressed for the ion trapped closer to the surface because of a larger density of states of surface Andreev bound states.
However, the observed mobility exhibits no depth dependence (Fig. \ref{3He-B_mobility}).
The lack of depth dependence is very puzzling and prevents a straightforward interpretation.
At the time when these experimental data were taken, Ikegami {\it et al.} interpreted that the lack of depth dependence is due to the zero transport cross section associated with the absence of density fluctuations arising from the Majorana condition \cite{Ikegami2013-2}.
However, this interpretation recently turned out to be incorrect;
Tsutsumi theoretically showed that surface Andreev bound states have some contributions to the mobility \cite{Tsutsumi2017}.
The calculated mobility are shown in Fig. \ref{3He-B_mobility} for depths of 21, 40, and 58 nm.
In the theory, the scattering process of surface Andreev bound states by the ions was calculated for the wavefunction given by Eq. (\ref{eq_SABS}) to obtain the mobility.
The temperature dependence quantitatively agrees with the experimental dependence.
Furthermore, the theoretical mobility does not have the depth dependence, which is also consistent with the experimental observation.
In fact, the theory shows that the mobility has almost no depth dependence up to $z \sim 2\kappa^{-1}$ \cite{Tsutsumi2017}.

As seen in Fig. \ref{3He-B_mobility}, the experimental mobility is suppressed from that in bulk $^3$He-B (blue line) at $T/T_c <$ 0.5.
This is due to the scattering of surface Andreev bound states, which becomes more dominant at lower temperatures.
Indeed, at $T/T_c=0.25$, more than half of the contribution to the mobility arises from surface Andreev bound states.

\begin{figure}
\begin{center}
\includegraphics[width=0.6\linewidth,keepaspectratio]{bound_states_bubble.eps}
\end{center}
\caption{\label{bound_states_bubble} (Color online)
(a) Energy spectra of bound states formed around a negative ion as a function of the quantum number of the angular momentum $l$. The energy levels of the bound states are linear in $l$. (The figure was provided by
courtesy of Y. Tsutsumi.)
(b) Transport cross section $\sigma_{\rm tr}$ of the negative ion for surface Andreev bound states at $z=0.5\xi$, $\xi$, and 1.5$\xi$, where $\xi$ is defined as $\xi=\kappa^{-1}/2$. Note that $\sigma_{\rm tr}$ becomes smaller as approaching the surface. Reprinted with permission from Ref. \cite{Tsutsumi2017}. Copyright (2013) by The American Physical
Society.}
\end{figure}

For the scattering of surface Andreev bound states, bound states that form around an ion play a crucial role \cite{Tsutsumi2017}.
When an ion is located in bulk $^3$He-B, bound states are formed around the ion with energies less than $\Delta_B$ [Fig. \ref{bound_states_bubble}(a)] \cite{Thuneberg1981-1,Thuneberg1981-2}.
These bound states are a microscopic realization of the surface Andreev bound states of $^3$He-B.
Because the bound states are formed in a small region around the ion, they have discrete energy levels, as shown in Fig. \ref{bound_states_bubble}(a).
The energy levels have linear dependence on the quantum number of the angular momentum $l$.
These bound states make a significant contribution to the scattering;
the surface Andreev bound states are resonantly scattered by the bound states around the ion when the energy of the surface Andreev bound state is close to one of the energy levels of the bound states formed around the ion.
In fact, the resonant scattering appears as sharp structures in the transport cross section $\sigma_{\rm tr}$ of the ion as a function of energy [Fig. \ref{bound_states_bubble}(b)].
Note that the resonant scattering is similar to that in $^3$He-A, where QPs are resonantly scattered by the bound states formed around the ion, which gives rise to the skew scattering of QPs \cite{Shevtsov2017,Shevtsov2016}.

The local density of states of surface Andreev bound states increases as approaching to the surface [$N_s(E,z) \sim e^{-2\kappa z}$].
However, the mobility does not exhibit depth dependence at $z \lesssim 2\kappa^{-1}$, as shown in Refs. \cite{Tsutsumi2017}.
The lack of depth dependence is caused by the reduction in $\sigma_{\rm tr}$ for surface Andreev bound states near the surface [see Fig. \ref{bound_states_bubble}(b)], which compensates for the increase in the density of states near the surface, making the mobility independent of depth.
The reduction in $\sigma_{\rm tr}$ arises from the strong mixing of the bound states around the ion and the surface Andreev bound states \cite{Tsutsumi2017}.

In Ref. \cite{Ikegami2013-2}, Ikegami {\it et al.} attributed the lack of depth dependence of mobility to the absence of scattering of the surface Andreev bound states from the ion, as a result of zero density fluctuations associated with the Majorana nature (i.e., $\left| u_{{\bf{k}}, \sigma } \right| = \left| v_{{\bf{k}}, \sigma } \right|$, where $u_{{\bf{k}}, \sigma }$ and $v_{{\bf{k}}, \sigma }$ are Bogoliubov coefficients).
However, the theory by Tsutsumi suggest that the surface Andreev bound states are scattered from the ion.
This can be explained as follows;
The hybridization of surface Andreev bound states and the bound states formed around the ion results in the violation of $\left| u_{{\bf{k}}, \sigma } \right| = \left| v_{{\bf{k}}, \sigma } \right|$, leading to nonzero transport cross section for the surface Andreev bound states.

If surface Andreev bound states are Majorana fermions, their dispersion relation should be gapless.
In the experiment described here, thermally excited surface Andreev bound states have an energy of $E/\Delta_B (\simeq T/1.76T_c) \sim $ 0.14 even at the lowest temperature ($T/T_c \sim 0.25$).
This means that the energy scale smaller than $E/\Delta_B \sim $ 0.14 cannot be resolved even at the lowest temperature, making it difficult to judge whether or not the dispersion relation is really gapless form the experiment.
Demonstrating the gapless dispersion by detecting surface Andreev bound states with lower energies is an important future challenge to verify the gapless dispersion relation.

\section{Conclusions and Future Prospects}

In this article, we reviewed that the mobility of the negative ion measured by Ikegami {\it et al.} \cite{Ikegami2013-2} enabled the direct experimental detection of surface Andreev bound states formed at a free surface.
In particular, the temperature and depth dependences of the experimental mobility have been quantitatively reproduced by a theoretical calculation that took into account the scattering of surface Andreev bound states formed at the specular surface \cite{Tsutsumi2017}.
The surface Andreev bound states formed at the specular surface should be Majorana fermions with a well-defined linear dispersion, and therefore the studies described in this article open the way for pursuing Majorana fermions under ideal conditions.

Proving the Majorana nature of surface Andreev bound states is an important next step.
One potential way of achieving this is to detect the Ising-like behavior of the spin of surface Andreev bound states \cite{Chung2009,Nagato2009}.
The Ising-like spin originates from the Majorana condition for the wavefunction of surface Andreev bound states, and may appear as anisotropic magnetic responses of the surface of $^3$He-B (and a $^3$He-B film) \cite{Chung2009,Nagato2009,Mizushima2012-1,Mizushima2012-2,Silaev2011} and the enhancement of the susceptibility at the surface \cite{Nagato2009,Mizushima2012-1,Mizushima2012-2}.
To detect the anisotropic magnetic response, experiments using ions trapped below a free surface have also been proposed \cite{Chung2009,Batulin2014}.

As mentioned above, detecting surface Andreev bound states with low energies is important for demonstrating a gapless linear dispersion.
Such low-energy states can in principle be detected by the mobility of ions at much lower temperatures as well as by observing power-law behaviors of thermodynamic quantities as a function of temperature, such as specific heat $C \propto T^2$ \cite{Mizushima2011} and superfluid mass fraction $\rho_s /\rho \propto 1-aT^3$ \cite{Wu2013}, for a $^3$He-B film at low enough temperatures ($a$ is a constant).
The detection of such low-energy states also allows us to study anisotropic gap-opening in a small magnetic field applied normal to the surface \cite{Chung2009,Mizushima2011} and the topological phase transition predicted in a magnetic field applied parallel to a surface \cite{Mizushima2012-1,Mizushima2012-2}.

An interesting approach to studying the properties of surface Andreev bound states in detail is to investigate nonlinear behaviors in the transport of an ion trapped below the surface.
When such an ion moves at a high velocity of order $\Delta_B/p_F$, the emission of surface Andreev bound states is expected by a mechanism in which the bound states formed around the ion escape into the Andreev surface bound states.
Although nonlinear properties have not been experimentally investigated for a negative ion,
 nonlinear transport has been measured for a positive ion \cite{Ikegami2014,Shiino2003}.
In the case of a positive ion, no depth dependence was found in the electric field--velocity relation \cite{Ikegami2014}, which has not yet been explained theoretically.
(The lack of depth dependence of the mobility observed for a positive ion \cite{Ikegami2013-2} has also not yet been explained.)
Experiments of nonlinear transport should be performed also for a negative ion.

The free surface of $^3$He-B is an exceptionally ideal playground for exploring Majorana fermions in condensed matter systems.
We hope that the Majorana nature will soon be demonstrated directly.
Unusual aspects of Majorana fermions might be uncovered by investigations at the free surface, which should be common to those expected in other condensed matter systems as well as in particle physics.

\section*{Acknowledgments}

We acknowledge Y. Tsutsumi, S.-B. Chung, T. Mizushima, O. Shevtsov, J. A. Sauls, and A. J. Leggett for illuminating discussions.
This work was partly supported by JSPS KAKENHI Grant Numbers JP24000007, JP26287084, and JP17H01145.


text:
\section{Introduction}\label{sec:introduction}
One of the distinctive characteristics of plasmas is that they have time and length scales that are typically well separated.
For example, in a 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 the 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 \citep{matthews94, waltz97, park99, cheng99, goswami05, brizard07, valentini07}. Such models could be ordered in a hierarchy that transitions all the way from kinetic to fluid theory \citep{dendy95, tronci14, tronci15, camporeale16}.
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 inertia-less 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 gyro-period are neglected) results into 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: $\omega\ll\Omega_i$, and that the wave spectrum is very anisotropic: $k_{||}\ll k_\perp$. 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 \emph{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 'gray 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 gray 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., \citep{howes11, valentini10, servidio12, franci15a}).
Indeed, it is now widely recognized that a kinetic treatment is necessary for studying solar wind turbulence, at small scales \citep{camporeale11, sahraoui12, boldyrev13, haynes14, vasconez14, haynes15}.
Therefore, it is legitimate to ask what is the range of validity of such models, and when 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 non-linear 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\rho_i=0.1$ to $k\rho_i=44$ (with $\rho_i$ the ion Larmor radius), and for a wide range of propagation angles $\theta$, ranging from $10^\circ$ to $85^\circ$; second, we present a robust and novel methodology to perform such 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 Section \ref{sec:methodology}, where we also define the errors that will be measured. Section \ref{sec:results} shows the comparison between hybrid, gyrokinetic, and fully-kinetic solutions for a wide range of cases. Conclusions are drawn in Section \ref{sec:conclusions}.

\section{Methodology and chosen parameters}\label{sec:methodology}
We study a plasma composed of electrons and protons, with mass ratio set to the physical value $m_i/m_e=1836$, in Cartesian coordinates $(\hat{i},\hat{j},\hat{k})$. As it is customary in linear theory, we study an uniform Maxwellian plasma with density $n_0$, in a constant magnetic field $\mathbf{B}=B_0\hat{k}$, that is perturbed with a small-amplitude plane wave of the form $\exp[i(\mathbf{k}\cdot \mathbf{x} -\omega t)]$, with wavevector $\mathbf{k}=(k_\perp,0,k_{||})$, and complex frequency $\omega=\omega_r + i\gamma$.\\
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 \citet{salem12}: $B_0 = 11$ nT, $T_e = 13.0$ eV (electron temperature), $T_p = 13.6$ eV (proton temperature), $n_0 = 9 $ cm$^{-3}$. The proton plasma beta is $\beta_p =0.4$, and the electron plasma to cyclotron frequency ratio is $\omega_{pe}/\Omega_e \simeq 88$. The ratio between Alfv\'{e}n and ion thermal velocities is $v_A/v_{ti}=1.57$, and the ratio between Alfv\'{e}n velocity and speed of light is $v_A/c=2.6\times10^{-4}$.\\
In all of the three cases considered (fully-kinetic, hybrid, and gyrokinetic), the linear wave solutions are obtained by finding the complex frequencies $\omega$ for which the determinant of a $3\times 3$ matrix $\mathbf{D}$ vanishes. $f(\omega,\mathbf{k})=\det(\mathbf{D})$ is a nonlinear function of both $\mathbf{k}$ and $\omega$, and the standard way to numerically solve the problem, for a fixed $\mathbf{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 $\omega$ varying $\mathbf{k}$. In other words, the solver can easily jump from a solution branch to another, whenever two solutions get close enough. The paper by \citet{podesta12} explains how to reduce the numerical error in calculating the elements of $\mathbf{D}$ in the case of the Vlasov-Maxwell system. \citet{podesta12} 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 analyze separately the real $(\Re)$ and imaginary $(\Im)$ parts of the function $f(\omega)$, and to identify the curves in the complex plane, where $\Re(f)=0$ and $\Im(f)=0$, that is where they change sign. Clearly, the intersections of such curves ($\Re(f)=\Im(f)=0$) are the roots $f(\omega)=0$, that can be used as initial guess for the Newton solver, to achieve a better accuracy. We have employed this method to identify \emph{all} the roots in a given region of the complex plane, in Section \ref{sec:def_err}\\
We refer to \citet{howes06} 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 $\varepsilon(\omega,\mathbf{k})$ is additive by species (see, e.g. \citet{stix62}):
\begin{equation}
\varepsilon(\omega,\bf{k}) = \bf{1} + \sum_s\chi_s(\omega,\bf{k}),
\end{equation}
where $\mathbf{1}$ is the unit dyadic, $\chi_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 $\chi_e$ derived within a two-fluid linear model. Here, we have employed the so-called warm plasma dispersion relation derived in \citet{swanson12}, 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 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).

\subsection{Diagnostics}\label{sec: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 \emph{oblique} propagation, for $\theta\leq 60^\circ$, and \emph{quasi-perpendicular} propagation, for $\theta>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 in considerations the overall spectral properties of the linear operator. Such approach is justified by the \emph{non-normality} of the Vlasov-Maxwell linear operator. We refer to, e.g., \citep{camporeale09, camporeale10, podesta10, ratushnaya14, friedman14} 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 \citet{trefethen05} and, for the case of the Vlasov-Maxwell system, in \citet{camporeale12}.
For instance, Figure \ref{fig:modes_80} shows a subset of the normal modes (obtained in the fully-kinetic model), in the complex plane, for $\theta=80^\circ$. The damping rate $\gamma$ and the real frequency $\omega_r$ are on horizontal and vertical axis, respectively (here and in all figures they are normalized to the ion cyclotron frequency $\Omega_i$). The different symbols are for $k\rho_i=0.1$ (blue dots), $k\rho_i=1$ (red circles), and $k\rho_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\'{e}n Wave (KAW) mode. This is particularly clear by inspection of Figure \ref{fig:modes_5_80}.
Here, we plot, in the top panel, the subset of normal modes that are confined in the complex plane within the domain $-2.5<\gamma<0$ and $0<\omega_r<20$.
The blue dots and red circles denote the fully-kinetic and the hybrid solutions, respectively. The Kinetic Alfv\'{e}n wave, that is the root of the dispersion relation that when tracked to low $k$ and low $\theta$ is connected to the Alfv\'{e}n wave mode, is denoted with a blue cross ($\omega_r = 0.47$, $\gamma=-1.29$). It is evident that, in this regime, KAW is surrounded by many other modes, some of them with similar frequency but lower damping rate. Nevertheless, KAW has been attributed a predominant role in some of the recent solar wind literature \citep{sahraoui09, salem12}. Even if one would like to adhere to the classical view where the least damped modes are more important than the heavily-damped one, it is clear that KAW is not special. The bottom panel of Figure \ref{fig:modes_5_80} shows the same modes, but now sorted as function of their damping rate per wave period $(|\gamma|/\omega_r)$ (horizontal axis), and their polarization $P$ (vertical axis). The polarization is here defined following \citet{hunana13}, as $P=\arg(E_y/E_x)/\pi$, 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, KAW is the most rapidly damped mode per wave period. A final remark on Figure \ref{fig:modes_5_80} is that one can see that, despite being in the electron range (i.e. $k\rho_i>1$), the hybrid model captures most of the normal modes quite accurately. It creates a spurious mode with $\omega_r=14$ and $\gamma=-0.08$ due to the inability to account for electron Landau damping.

\subsubsection{Definition of errors}\label{sec:def_err}
We define the following two measures of error, that assess the deviation of the frequency $\omega$ (obtained either with hybrid or gyrokinetic), with respect to the correct value obtained with the Vlasov-Maxwell theory, denoted as $\omega_{VM}$:
\begin{equation}
 \varepsilon_\omega = \frac{|\omega_{VM}-\omega|}{|\omega_{VM}|}\label{err_1}
\end{equation}
\begin{equation}
 \varepsilon_{wave} = ||\Re\left({\exp[-i\omega_{VM} t]}\right) - \Re\left({\exp[-i\omega t]}\right) ||_2\label{err_2}
\end{equation}
The error $\varepsilon_\omega$ in Eq.(\ref{err_1}) 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 $\varepsilon_{wave}$ in Eq.(\ref{err_2}), we define a time interval $t=[0, 2\pi/\Re(\omega_{VM})]$, that covers one wave period (with respect to $\omega_{VM}$), and 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 $\omega$ will differ, during one wave period, from the reference wave with frequency $\omega_{VM}$, by modeling its actual time evolution.\\
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 $\varepsilon$-pseudospectrum $\Lambda_\varepsilon(\mathbf{D})$, for a fixed $\mathbf{k}$:
\begin{equation}
 \Lambda_\varepsilon(\mathbf{D})=\{z\in\mathbb {C}: |\det(\boldsymbol{D}(z))| \leq \varepsilon \}.
\end{equation}
Clearly, the pseudospectrum reduces to the standard spectrum (set of eigenvalues) for $\varepsilon\rightarrow 0$.
For instance, Figures \ref{fig:pseudo} and \ref{fig:pseudo_k_1_85} compare the fully-kinetic (left), hybrid (middle), and gyrokinetic (right) pseudospectra for the case $k\rho_i=1$ and $\theta=80^\circ$ (Figure \ref{fig:pseudo}) or $\theta=85^\circ$ (Figure \ref{fig:pseudo_k_1_85}), in the region $-0.5<\gamma<0$, and $0<\omega_r<1$. The color plot is in logarithmic scale and denotes the value of $\det(\mathbf{D})$.
One can visually identify the eigenvalues as the points where the concentric curves are centered.
In order to measure quantitatively the error in the pseudospectra, we proceed as follows.
We compute, for a fixed value of $k=|\mathbf{k}|$ the function $f(\omega)=\det(\mathbf{D(\omega)})$ on a $500\times 500$ grid covering the domain
$0\leq\omega_r\leq5k$, and $-0.5k\leq\gamma\leq0$. Because $f$ can vary 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(\omega_{VM})$ is small, disregarding regions where $f$ is large. Therefore, we map $f(\omega)$ into a new function $\phi(\omega)$ through a sigmoid:
\begin{equation}
 \phi(\omega) = 1+ \frac{9f}{1+f}.
\end{equation}
The new function $\phi(\omega)$ is defined between 1 (for $f=0$) and 10 (for $f\rightarrow\infty$).
Finally, we define the error $\varepsilon_D$ as:
\begin{equation}
 \varepsilon_D = \left\|\frac{ (\phi(\omega) - \phi(\omega_{VM}))}{\phi(\omega_{VM})}\right\|_2,\label{eps_D}
\end{equation}
where the L-2 norm is normalized on the number of grid points.

\subsection{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_\perp$. 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 \ref{fig:validity_GK} shows the range of validity of the gyrokinetic model. The symbols denote the angle of propagation $\theta$ for a given value of perpendicular wavevector $k_\perp\rho_i$, that results in a relative error on the Kinetic Alfv\'{e}n Wave real frequency of $10\%$ (squares) and $30\%$ (diamonds). One can appreciate that for $k_\perp\rho_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 \ref{fig:validity_GK_2}, that shows the dispersion relation for the Alfv\'{e}n/ion-cyclotron wave. When plotted as $\omega_r/k_{||}v_A$ as a function of $k_\perp\rho_i$, the dispersion relation becomes independent of angle, in gyrokinetics (shown in blue). The red, yellow, and purple lines denote the fully-kinetic solution for angles of propagation $\theta = 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 \citep{howes11}, on the basis that some observational data suggest that the propagation angle in the solar wind averages at about 88 degrees, and across the examined spectra never fall below 60 degrees \citep{sahraoui10, narita10}.
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 \citep{osman09, kiyani13, lion16}. Furthermore, the development of wave vector 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.

\section{Results}\label{sec:results}
Historically, different communities in plasma physics have adopted different nomenclatures for linear waves. Here we follow the nomenclature adopted by \citet{gary_book}. We will refer to Alfv\'{e}n/ion-cyclotron 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 pseudo-spectra in the quasi-perpendicular regime.

\subsection{Alfv\'{e}n/ion-cyclotron (AIC)}
Figure \ref{fig:disp_rel_AIC} shows the dispersion relation for the AIC mode, for angles of propagation $\theta = 10^\circ, 30^\circ, 60^\circ$. Blue, red and green curves are for fully-kinetic, hybrid and gyrokinetic, respectively (the same notation will be used for all following figures). The real frequency $\omega_r$ is shown in the top panel, and the damping rate in the bottom panel. In both panels $k_\perp\rho_i$ is on the horizontal axis (logarithmic scale is used on all axis). The three different propagation angles can be identified from the starting values of $k_\perp\rho_i$. Since all the curves are plotted from $k\rho_i=0.1$, the initial perpendicular wavevector increases with increasing $\theta$. For $\theta=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 $\theta=60^\circ$ solution is that a branch bifurcation happens around $k_\perp\rho_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_\perp\rho_i<1$. The gyrokinetic solutions are generally less damped than the correct fully-kinetic solutions (bottom panel).
The errors for the AIC mode solutions are shown in Figures \ref{fig:error_AIC_10}, \ref{fig:error_AIC_30}, \ref{fig:error_AIC_60}, for the cases $\theta=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 $\varepsilon_\omega$ and $\varepsilon_{wave}$, respectively. One can notice that the errors for the hybrid solutions are lower than the corresponding errors for the gyrokinetic model, for $\theta=10^\circ$ and $30^\circ$. However, which model performs better for $\theta=60^\circ$ depends on $k_\perp\rho_i$
The fact that the gyrokinetic errors in Figure \ref{fig:error_AIC_60} for $\theta=60^\circ$ do not grow very large, and almost saturate for $k_\perp\rho_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 \ref{fig:disp_rel_AIC_perp} shows the dispersion relation for AIC for $\theta = 80^\circ, 85^\circ, 89^\circ$ (from left to right). Interesting, even at such large propagation angle, the AIC resonate 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 (bottom panel). Moreover, in the $\theta = 89^\circ$ case, the saturated value of the hybrid solution is not correct (being larger than $\omega/\Omega_i =1 $). The respective errors are plotted in Figures \ref{fig:error_AIC_80}, \ref{fig:error_AIC_85}, and \ref{fig:error_AIC_89}. 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).

\subsection{Ion-acoustic (IA)}
The dispersion relation for the IA mode is shown in Figure \ref{fig:disp_rel_IA}, with the same format as in the previous Section. One can distinguish the following features. Regarding the real frequency (top panel), there is again a good agreement between hybrid and fully-kinetic; the gyro-kinetic frequency also presents a similar agreement, but only for small $k$ values. Damping rates present larger mismatches. Notably, the dispersion relation for the case $\theta=60^\circ$ presents an 'elbow' for $k_\perp\rho_i\sim 1$ that is captured by the hybrid model.
We show the errors $\varepsilon_\omega$ and  $\varepsilon_{wave}$ for $\theta=10^\circ$, $\theta=30^\circ$, $\theta=60^\circ$ respectively in Figures \ref{fig:error_IA_10}, \ref{fig:error_IA_30}, and \ref{fig:error_IA_60} (same format as previous Section).
The errors are now generally smaller for gyrokinetic at small $k_\perp\rho_i$, and increasing with wavector, and almost constant in value fo rhte hybrid solutions.
The dispersion relation plots for quasi-perpendicular angles are shown in Figure \ref{fig:disp_rel_IA_perp}. 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 $\theta=89^\circ$.
The corresponding errors are shown in Figures \ref{fig:error_IA_80}, \ref{fig:error_IA_85}, and \ref{fig:error_IA_89}. It is interesting to notice how the relative accuracy of the two models strongly depend on the wavevector. For instance, at $\theta = 80^\circ$ the gyrokinetic solution is about two order of magnitudes more accurate than the hybrid for $k_\perp\rho_i = 0.1$, but starts to become less accurate for
$k_\perp\rho_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.

\subsection{Fast mode}
It is well-known that the fast mode solution is factored out from the gyrokinetic model. In this sense, a comparison of gyrokinetics with the other the two models is meaningless.
In Figure \ref{fig:disp_rel_FW} we compare the dispersion relations of the fully-kinetic, the hybrid, and the two-fluids models for $\theta=10^\circ$, $\theta=30^\circ$, $\theta=60^\circ$. The dispersion relation for the latter follows the derivation of \citet{zhao14}. The three solutions are in a very good agreement for the real frequency (top panel), where a slight deviation with respect to the correct solution occurs only at large $k_\perp$. Interestingly, the two-fluids solution is a very good proxy for both hybrid and full-kinetic solutions, in this regime.
The bottom panel shows the corresponding damping rate (the two-fluids 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\rho_i\leq 1$, where the red and blue curves have a good agreement.

\subsection{Quasi-perpendicular propagation}
In this section we perform the analysis of the error $\varepsilon_D$ of Eq. (\ref{eps_D}). 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 \ref{fig:error_D} shows the error $\varepsilon_D$ for angles of propagation $\theta=60^\circ$, $\theta=70^\circ$, $\theta=80^\circ$, and $\theta=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_\perp\rho_i=1$, while the hybrid model error tend to increase for $k_\perp\rho_i<1$ and to decrease for $k_\perp\rho_i>1$. Overall, the hybrid model performs much better than the gyrokinetic. This is not surprising because, as it was already evident by inspection of Figures \ref{fig:pseudo} and \ref{fig:pseudo_k_1_85}, 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 $\varepsilon_D$ is defined and constructed in a rather arbitrary way. Although this is certainly true, we think that it is still very informative on how much a reduced model is able to reproduce the correct linear evolution of the system. as we mentioned, the ability of gyrokinetics to describe solar wind turbulence in the regime $k_\perp\gg k_{||}$ has been stressed out in many recent papers \citep{schekochihin09, tenbarge13}.
Although this might certainly be true, we believe that it is informative to know what is left out from the model, in terms of normal modes. In Figure \ref{fig:modes_contours} we compare the three models (fully-kinetic to the left, hybrid in the middle, gyrokinetic to the right) for the case $k\rho_i=1$, and $\theta=85^\circ$. Here, we do not show the normal modes, but the curves $\Re(f)=0$ (black), and $\Im(f)=0$ (red) in the domain $-0.2<\gamma<0$ and $0<\omega_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: gyrokinetics 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).

\section{Conclusions}\label{sec: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 both dispersion relations of single normal modes, for angles of propagation $\theta = 10^\circ, 30^\circ, 60^\circ, 80^\circ, 85^\circ, 89^\circ$ and pseudospectra, for   $\theta = 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 that the gyrokinetic model has been employed also outside its nominal range of validity ($k_{||}\ll k_\perp$).
We note that the Vlasov-Maxwell linear theory has recently been investigated by \citet{sahraoui12}, 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 $\theta=89^\circ$. Interestingly, the errors decrease in the region when the waves resonate with ions, and the dispersion curves saturate. One major mismatch of the hybrid model happens for $\theta = 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_\perp\rho_i>1$ However, also for quasi-perpendicular propagation, the errors increase abruptly as soon as the waves become resonant with ions.
In conclusion, it appears that the extension of the gyrokinetic model to the electron regime ($k_\perp\rho_i>1$) is not justified, for angles of propagation smaller than about $85^\circ$. This is also apparent from looking at the analysis of the error based on the pseudospectra. Indeed, it appears that the hybrid model should be preferred, for its ability to represent most of the normal modes (certainly in the ion regime $k_\perp\rho_i<1$ and, to a certain extent, in the electron regime as well). 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, still retaining the fully-kinetic nature of the solution (see, e.g. \citep{degond12, chen14, tronci15}).

\newpage

\begin{figure}
  \includegraphics[width=10cm]{modes_80.eps}
  \caption{Normal modes derived with the fully-kinetic model. Damping rate ($\gamma$) and real frequency $(\omega_r)$ are on horizontal and vertical axis, respectively. The angle of propagation is $\theta=80^\circ$. Blue dots, red circles, and black square denote $k\rho_i=0.1$, 1, 5, respectively.}
\label{fig:modes_80}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{modes_5_80.eps}
  \caption{Normal modes for $\theta=80^\circ$ and $k\rho_i=5$ for fully-kinetic (blue dots) and hybrid (red circles).
  Top panel: solutions in the complex plane $(\gamma,\omega_r)$. Bottom panel: solutions plotted as functions of $|\gamma|/\omega_r$ (horizontal axis) and polarization P (vertical axis). In both panels the Kinetic Alfv\'{e}n Wave (KAW) is denoted with a blue cross.}
\label{fig:modes_5_80}
\end{figure}

\begin{figure}
  \includegraphics[width=15cm]{pseudo.eps}
  \caption{Pseudospectra in the complex plane $(\gamma,\omega_r)$ for $k\rho_i = 1$, and $\theta = 80^\circ$. Left: full-kinetic; middle: hybrid; right: gyrokinetic. The color plot represents $f(\omega)=\det(\bf{D})$ in logarithmic scale.}
\label{fig:pseudo}
\end{figure}

\begin{figure}
  \includegraphics[width=15cm]{pseudo_k_1_85.eps}
  \caption{Pseudospectra in the complex plane $(\gamma,\omega_r)$ for $k\rho_i = 1$, and $\theta = 85^\circ$. Left: full-kinetic; middle: hybrid; right: gyrokinetic. The color plot represents $f(\omega)=\det(\bf{D})$ in logarithmic scale.}
\label{fig:pseudo_k_1_85}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{validity_GK.eps}
  \caption{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_\perp \rho_i$ and the angle of propagation $\theta$ are on horizontal and vertical axis. The model is valid above the symbols. }
\label{fig:validity_GK}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{validity_GK_2.eps}
  \caption{Range of validity of the gyrokinetic model. Dispersion relation for Alfv\'{e}n/ion-cyclotron wave, plotted as $\omega_r/k_{||}v_A$ as a function of $k_\perp\rho_i$. The blue line denotes the gyrokinetics solution, which is independent of the angle of propagation. The red, yellow, and purple lines denote the fully-kinetic solution for angles of propagation $\theta = 80^\circ,85^\circ,89^\circ$. }
\label{fig:validity_GK_2}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{disp_rel_AIC.eps}
  \caption{Dispersion relation for Alfv\'{e}n/ion-cyclotron waves for angles of propagation $\theta=10^\circ, 30^\circ, 60^\circ$ (curves from left to right). Blue, red, and green lines are for fully-kinetic, hybrid, and gyrokinetic, respectively. Top panel: real frequency $\omega_r$; bottom panel: damping rate $\gamma$. $k_\perp\rho_i$ is on the horizontal axis. All axis are in logarithmic scale.}
\label{fig:disp_rel_AIC}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_10.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=10^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_10}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_30.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=30^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_30}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_60.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=60^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_60}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{disp_rel_AIC_perp.eps}
  \caption{Dispersion relation for Alfv\'{e}n/ion-cyclotron waves for angles of propagation $\theta=80^\circ, 85^\circ, 89^\circ$ (curves from left to right). Blue, red, and green lines are for fully-kinetic, hybrid, and gyrokinetic, respectively. Top panel: real frequency $\omega_r$; bottom panel: damping rate $\gamma$. $k_\perp\rho_i$ is on the horizontal axis. All axis are in logarithmic scale.}
\label{fig:disp_rel_AIC_perp}
\end{figure}
\clearpage

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_80.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=80^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_80}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_85.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=85^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_85}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_AIC_89.eps}}
  \caption{Errors for Alfv\'{e}n/ion-cyclotron wave with $\theta=89^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_AIC_89}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{disp_rel_IA.eps}
  \caption{Dispersion relation for ion-acoustic waves for angles of propagation $\theta=10^\circ, 30^\circ, 60^\circ$ (curves from left to right). Blue, red, and green lines are for fully-kinetic, hybrid, and gyrokinetic, respectively. Top panel: real frequency $\omega_r$; bottom panel: damping rate $\gamma$. $k_\perp\rho_i$ is on the horizontal axis. All axis are in logarithmic scale.}
\label{fig:disp_rel_IA}
\end{figure}

\clearpage

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_10.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=10^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_10}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_30.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=30^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_30}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_60.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=60^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_60}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{disp_rel_IA_perp.eps}
  \caption{Dispersion relation for ion-acoustic waves for angles of propagation $\theta=80^\circ, 85^\circ, 89^\circ$ (curves from left to right). Blue, red, and green lines are for fully-kinetic, hybrid, and gyrokinetic, respectively. Top panel: real frequency $\omega_r$; bottom panel: damping rate $\gamma$. $k_\perp\rho_i$ is on the horizontal axis. All axis are in logarithmic scale.}
\label{fig:disp_rel_IA_perp}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_80.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=80^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_80}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_85.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=85^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_85}
\end{figure}

\begin{figure}
  \centerline{\includegraphics[width=10cm]{error_IA_89.eps}}
  \caption{Errors for ion-acoustic wave with $\theta=89^\circ$, as functions of $k_\perp\rho_i$. Red and green curves are for hybrid and gyrokinetics, respectively. Solid lines denote $\varepsilon_\omega$, and dashed lines denote $\varepsilon_{wave}$.}
\label{fig:error_IA_89}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{disp_rel_FW.eps}
  \caption{Dispersion relation for fast waves for angles of propagation $\theta=10^\circ, 30^\circ, 60^\circ$ (curves from left to right). Blue, red, and green lines are for fully-kinetic, hybrid, and two-fluid, respectively. Top panel: real frequency $\omega_r$ in logarithmic scale; bottom panel: damping rate $\gamma$ in logarithmic scale. $k_\perp\rho_i$ is on the horizontal axis in logarithmic scale.}
\label{fig:disp_rel_FW}
\end{figure}

\begin{figure}
  \includegraphics[width=10cm]{error_D.eps}
  \caption{Error $\varepsilon_D$, as a function of $k_\perp\rho_i$, for angles of propagation $\theta=60^\circ, 70^\circ, 80^\circ, 85^\circ$ (blue, red, yellow, and magenta). Dashed and solid lines are for gyrokinetics and hybrid, respectively.}
\label{fig:error_D}
\end{figure}

\begin{figure}
  \includegraphics[width=15cm]{modes_contours.eps}
  \caption{Real (black) and imaginary (red) part of $f(\omega)=\det(\bf{D})$ for $\theta=85^\circ$, $k\rho_i=1$. Left panel: fully-kinetic; middle: hybrid; right: gyrokinetic. The intersection between black and red curves identify normal modes.}
\label{fig:modes_contours}
\end{figure}

\clearpage

\bibliographystyle{jpp}

\newpage
\clearpage

redpajama_set_name: RedPajamaArXiv