# Structure and evolution of magnetohydrodynamic solitary waves with Hall and finite Larmor radius effects

E. Bello-Benítez, G. Sánchez-Arriaga, T. Passot, D. Laveder, E. Siminos

► To cite this version:

E. Bello-Benítez, G. Sánchez-Arriaga, T. Passot, D. Laveder, E. Siminos. Structure and evolution of magnetohydrodynamic solitary waves with Hall and finite Larmor radius effects. Physical Review E, 2019, 99 (2), 10.1103/PhysRevE.99.023202. [https://arxiv.org/abs/1902.02320]

HAL Id: hal-02303360

https://hal.science/hal-02303360

Submitted on 10 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Structure and evolution of magnetohydrodynamic solitary waves with Hall and finite Larmor radius effects

E. Bello-Benítez,¹,* G. Sánchez-Arriaga,¹,† T. Passot,²,‡ D. Laveder,²,§ and E. Siminos³,||

¹Bioengineering and Aerospace Engineering Department, Universidad Carlos III de Madrid, Leganés, Spain

²Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France

³Department of Physics, University of Gothenburg, 412 96 Göteborg, Sweden

<figure><img src="image_2.png" /></figure>

(Received 22 October 2018; published 8 February 2019)

Nonlinear and low-frequency solitary waves are investigated in the framework of the one-dimensional Hall-magnetohydrodynamic model with finite Larmor effects and two different closure models for the pressures. For a double adiabatic pressure model, the organization of these localized structures in terms of the propagation angle with respect to the ambient magnetic field $\theta$ and the propagation velocity $C$ is discussed. There are three types of regions in the $\theta-C$ plane that correspond to domains where either solitary waves cannot exist, are organized in branches, or have a continuous spectrum. A numerical method valid for the two latter cases, which rigorously proves the existence of the waves, is presented and used to locate many waves, including bright and dark structures. Some of them belong to parametric domains where solitary waves were not found in previous works. The stability of the structures has been investigated by performing a linear analysis of the background plasma state and by means of numerical simulations. They show that the cores of some waves can be robust, but, for the parameters considered in the analysis, the tails are unstable. The substitution of the double adiabatic model by evolution equations for the plasma pressures appears to suppress the instability in some cases and to allow the propagation of the solitary waves during long times.

DOI: 10.1103/PhysRevE.99.023202

# I. INTRODUCTION

Solitary waves of various types are commonly observed in collisionless heliospheric plasmas. A convincing observational evidence of large amplitude electromagnetic solitary wave propagating in the terrestrial environment was provided by Cluster multisatellite data near the magnetopause boundary [1]. The soliton, whose size is a few inertial lengths, is of slow type and is relatively stable as it displayed very similar shapes when observed from two satellites at two different physical locations. Other types of nonlinear waves in the form of fast magnetosonic shocklets are also observed with the Cluster satellites near the earth bow shock [2]. Compressive solitary structures or shocks are identified even further in the slow solar wind [3]. Various structures in the form of single nonlinear Alfvén wave cycles, discontinuities, magnetic decreases, and shocks embedded in the turbulence of high-speed solar wind streams are reviewed in Ref. [4]. Magnetic humps or holes in total pressure balance, either in the form of isolated structures or in wave trains, are also commonly observed in planetary magnetosheath or in the solar wind [5]. They are often attributed to the saturation of the mirror instability. The latter being subcritical [6], isolated magnetic holes can also

be observed below the threshold of the mirror instability. Magnetic humps, on the other hand, often require sufficiently large temperature anisotropy. These structures are clearly different from slow or fast modes as they are nonpropagating in the plasma rest frame. Their propagation velocity is difficult to measure precisely, however, so that some uncertainty subsists in their identification. A complete determination of the various hydrodynamic as well as electromagnetic fields could permit one to alleviate the ambiguity, but this also remains a difficult observational task. Even though these nonlinear structures are observed in almost collisionless plasmas at scales of the order of a few ion Larmor radii, fluid modeling, possibly accounting for ion finite Larmor radius effects, appears to be sufficient to reproduce their main properties. Their amplitude is large, however, and an important challenge is to describe them as solutions of the fully nonlinear extended fluid equations rather than small-amplitude asymptotic models.

Theoretical works on one-dimensional, localized, and traveling structures have contributed to the understanding of the propagation of nonlinear and low-frequency waves in plasmas. In the small-amplitude limit, these waves are governed by standard integrable equations such as the Korteweg–de Vries (KdV) [7], the modified KdV (MKDV) [8], the derivative nonlinear Schrödinger equation (DNLS) [9–11], and the triple-degenerate DNLS equation [12]. Some of them admit solitonic solutions, and relations between their propagation velocity and their amplitudes exist (see Ref. [13] and references therein). These small-amplitude asymptotic equations are also well suited to address questions related to perturbations of these solitary structures, such as the nontrivial effect

<footnote>*ebello@ing.uc3m.es</footnote>

<footnote>†gonzalo.sanchez@uc3m.es</footnote>

<footnote>‡thierry.passot@oca.eu</footnote>

<footnote>§dimitri.laveder@oca.eu</footnote>

<footnote>||evangelos.siminos@physics.gu.se</footnote>
of dissipation on Alfvén solitons [14], or to the role played by non-Maxwellian distribution functions on the shape and existence of solitons [15–17]. However, for finite amplitude, these localized structures should be studied in a more general framework, such as the magnetohydrodynamic models extended to include dispersive and/or dissipative effects. After assuming the traveling wave ansatz, the system of partial differential equations becomes a set of ordinary differential equations that can be used to investigate the existence of solitary waves and discontinuities. This technique has been used to study the structure of intermediate shock waves in the resistive-magnetohydrodynamics (MHD) [18], the resistive Hall-MHD [19], and the Hall-MHD with a double-adiabatic pressure tensor [20] systems, as well as rotational discontinuities in the Hall-MHD model with finite-Larmor-radius (FLR) and scalar pressure [21]. Exact solitary wave solutions in the Hall-MHD model for cold [22] and warm plasmas with scalar [23] and double-adiabatic pressure models [13] have been also found.

In the case of the Hall-MHD model with a double adiabatic pressure tensor, the traveling wave ansatz leads to a pair of coupled ordinary differential equations that governs the normalized components of the magnetic field normal to the propagation direction, called $b_y$ and $b_z$. Such a system has a hamiltonian structure and is reversible, i.e., solutions are invariant under the transformation $(\zeta, b_y, b_z) \to (-\zeta, -b_y, b_z)$, with $\zeta$ the independent variable. Adding FLR effects does not change the reversible character of the dynamical system, but it increases the effective dimension from two to four [24]. A numerical study of the existence of solitary waves in the parametric domain where the upstream state is a saddle-center was presented in Ref. [24]. Some geometrical arguments based on properties of the dynamical system, including its reversible character and the stability of the upstream state, were used in Ref. [20] to anticipate the organization of the solitary waves in parameter space by using well-known results for homoclinic orbits in reversible systems [25]. However, closing the equations with a double-adiabatic static model for the pressures relies on a set of assumptions, like neglecting the heat fluxes and using an ideal Ohms’ law. Retaining the Hall term in the Faraday equation then implies that the conservation of the energy is not guaranteed and, even for isotropic upstream conditions, that the background plasma state can be unstable. Both issues are overcome with and without FLR effects by using appropriate evolution equations for the pressures [26] (still assuming zero heat fluxes). Performing the analysis on the resulting system, however, becomes a much more difficult task.

Due to these circumstances, we split this work in two different phases. First, by taking advantage of the lower dimension of the phase space, we found solitary wave solutions in the FLR-Hall-MHD model with a double adiabatic pressure tensor. This phase covers the presentation of the model and its main properties in Sec. II and a numerical scheme to prove the existence of solitary waves rigorously in Sec. III. Several types of orbits and their main physical features are also presented. In a second phase, which is presented in Sec. IV, the stability of these waves is studied by using two different types of numerical simulations of the FLR-Hall-MHD model. It is shown that, if the equations are closed with the double

adiabatic model, some solitary waves have a highly unstable core, and others can propagate longer times but eventually are destroyed due to the instability of the background plasma state. On the other hand, simulations closed with evolution equations for the pressures shows that some solutions found in the first phase are very robust and their cores seem to propagate without significant deformation indefinitely. This result justifies our approach and highlights the existence of stable localized structures in the FLR-Hall-MHD model.

## II. THE FLR-HALL MHD MODEL

This section presents the dynamical system that governs the solitary waves in the FLR-Hall-MHD model closed by a double adiabatic pressure tensor. Although it was already found in a previous work [24], the steps in the derivation are summarized here in order to correct some errors and discrepancies. Cumbersome calculations are moved to Appendix A. The analysis is carried out in the framework of the Hall-MHD model with finite-Larmor-radius (FLR) effects. Mass density $\rho$, plasma (i.e., ion) flow velocity $v$, and magnetic field $B$ are governed by

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \quad (1) $$

$$ \frac{\partial}{\partial t}(\rho \mathbf{v}) + \nabla \cdot \left[ \rho \mathbf{v}\mathbf{v} + \bar{\mathbf{P}}_i + p_e \bar{\mathbf{I}} + \frac{1}{4\pi} \left( \frac{1}{2} B^2 \bar{\mathbf{I}} - \mathbf{B}\mathbf{B} \right) \right] = 0, \quad (2) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left[ \mathbf{v} \times \mathbf{B} - \frac{m_i c}{4\pi e \rho} (\nabla \times \mathbf{B}) \times \mathbf{B} \right], \quad (3) $$

where $m_i$ is the ion mass, $c$ the speed of light, $e$ the electron charge, and $\bar{\mathbf{I}}$ the identity tensor. We assumed that the electron pressure is isotropic and follows an isothermal equation of state $p_e = \rho v_{se}^2$ with $v_{se}^2 = k_B T_e / m_i$ the electron contribution to the ion-acoustic velocity, $T_e$ the electron temperature, and $k_B$ Boltzmann’s constant. The ion pressure tensor $\bar{\mathbf{P}}_i$ is written as

$$ \bar{\mathbf{P}}_i = \bar{\mathbf{P}}_{i,0}^{(0)} + \bar{\mathbf{P}}_{i,1}^{(1)} + \bar{\mathbf{P}}_{i,2}^{(1)} + \bar{\mathbf{P}}_{i,3}^{(1)}, \quad (4) $$

where the tensor with superscript 0 represents the gyrotropic contribution and reads

$$ \bar{\mathbf{P}}_{i}^{(0)} = p_{\parallel} \mathbf{e}_{b} \mathbf{e}_{b} + p_{\perp} (\bar{\mathbf{I}} - \mathbf{e}_{b} \mathbf{e}_{b}) \equiv \bar{\mathbf{P}}_{i,\parallel}^{(0)} + \bar{\mathbf{P}}_{i,\perp}^{(0)}, \quad (5) $$

with $p_{\parallel}$ and $p_{\perp}$ the parallel and perpendicular pressures and $\mathbf{e}_b = \mathbf{B}/B$ the unit vector along the magnetic field. Tensors with superscript 1 in Eq. (4) represent the FLR corrections and are given by [27,28]

$$ \bar{\mathbf{P}}_{i,1}^{(1)} = \frac{1}{\Omega_{ci}} \left[ \frac{1}{4} \mathbf{e}_b \times (\nabla \mathbf{v} + \nabla \mathbf{v}^T) \cdot \bar{\mathbf{P}}_{i,\perp}^{(0)} + \text{transp.} \right], \quad (6) $$

$$ \bar{\mathbf{P}}_{i,2}^{(1)} = -\frac{1}{\Omega_{ci}} [\mathbf{e}_b(\nabla \times \mathbf{v}) \cdot \bar{\mathbf{P}}_{i,\perp}^{(0)} + \text{transp.}], \quad (7) $$

$$ \bar{\mathbf{P}}_{i,3}^{(1)} = \frac{2}{\Omega_{ci}} [\mathbf{e}_b(\bar{\mathbf{P}}_{i,\parallel}^{(0)} \cdot \nabla) \times \mathbf{v} + \text{transp.}], \quad (8) $$

where $\Omega_{ci} = eB/(m_i c)$ is the local ion gyro frequency and the notation “+transp” means that one should sum the transpose of the tensor immediately to the left in the square bracket. The
equations are completed with the following double-adiabatic model for the equations of state

$$ \frac{p_{\parallel} B^2}{\rho^3} = \text{const}, \qquad (9) $$

$$ \frac{p_{\perp}}{\rho B} = \text{const.} \qquad (10) $$

Hereafter, subscript 0 will be used to denote the unperturbed variables. Therefore, $\rho_0$, $B_0$, $p_{\parallel 0}$, and $p_{\perp 0}$ correspond to the values of $\rho$, $B$, $p_{\parallel}$, and $p_{\perp}$ upstream from the solitary wave. We also now introduce a Cartesian frame of reference with the x axis along the propagation direction of the wave, and the y and z axes chosen such that the upstream magnetic field has no component in the y direction. Such a frame is linked to the solitary wave and moves at velocity $v_{x0}$ with respect to the unperturbed plasma. In the upstream region, i.e., at $x \to +\infty$, plasma velocity and magnetic field then read

$$ \mathbf{v}(x \to +\infty) = v_{x0} \mathbf{e}_x, \qquad (11) $$

$$ \mathbf{B}_0 = B_0(\cos \theta \mathbf{e}_x + \sin \theta \mathbf{e}_z), \qquad (12) $$

with $\mathbf{e}_x$, $\mathbf{e}_y$, and $\mathbf{e}_z$ unit vectors along the axes of the Cartesian frame. Therefore, the solitary wave propagates along the positive (negative) x direction for $v_{x0} < 0$ ($v_{x0} > 0$). We will consider the case $v_{x0} < 0$ and will use the wave velocity $C = -v_{x0}$.

If the analysis is restricted to stationary ($\partial/\partial t = 0$) and one-dimensional waves ($\partial/\partial y=\partial/\partial z=0$), then one finds that $B_x$ is constant ($B_x = B_0 \cos \theta$) and the FLR-Hall-MHD model becomes the set of ordinary differential equations

$$ \frac{d\xi}{d\hat{x}} = \mathbf{f}(\xi). \qquad (13) $$

The state vector of this dynamical system is five-dimensional, $\xi = [u_x \ u_y \ u_z \ b_y \ b_z]^T$, and it involves the normalized velocity $\mathbf{u} = u_x \mathbf{e}_x + u_y \mathbf{e}_y + u_z \mathbf{e}_z$ and magnetic field components $\mathbf{b} = b_y \mathbf{e}_y + b_z \mathbf{e}_z$ with $\mathbf{u} = \mathbf{v}/v_{x0}$ and $b_{y,z} = B_{y,z}/B_0 \sin \theta$. The explicit form of the vector flow $\mathbf{f} = [f_u \ f_b]$ in Eq. (13) can be found in Appendix A with $f_b$ and $f_u$ given by Eqs. (A5) and (A32).

The independent variable in Eq. (13) is the normalized length $\hat{x} = x/\ell$, with

$$ \ell = \frac{v_A^2 \cos \theta}{\Omega_{ci,0} v_{x0}}, \qquad (14) $$

$v_A = \sqrt{B_0^2/4\pi \rho_0}$ the Alfvén velocity, and $\Omega_{ci,0} = eB_0/(m_i c)$ the upstream ion cyclotron frequency. In Ref. [24] a $\cos \theta$ factor was missed. The dynamical system involves five parameters:

$$ \theta, \quad M_A = \frac{v_A^2}{v_{x0}^2}, \quad M_e = \frac{v_{se}^2}{v_{x0}^2}, \quad M_i = \frac{v_{\perp}^2}{v_{x0}^2}, \quad a_p = \frac{p_{\parallel 0}}{p_{\perp 0}}, \qquad (15) $$

$\theta$ being the angle between the propagation direction and the ambient unperturbed magnetic field and $v_{\perp}^2 = p_{\perp 0}/\rho_0$. This work investigates the effect of $C/V_A$ and $\theta$ on the properties of the solitary waves and will fix the other parameters according to the two cases shown in Table I.

<table><caption>TABLE I. Solitary wave parameters.</caption><thead><tr><th>Case</th><th>a<sub>p</sub></th><th>v<sub>&perp;</sub>/v<sub>A</sub></th><th>v<sub>se</sub>/v<sub>&perp;</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0.4</td><td>1</td></tr><tr><td>2</td><td>1.5</td><td>1.2</td><td>0.3</td></tr></tbody></table>

### A. Properties of the FLR-Hall MHD dynamical system

Before discussing interesting physical features of the solitary waves in Sec. III, we now summarize some purely mathematical results that are essential in order to organize the waves in the parameter space and design numerical algorithms to compute them and prove their existence rigorously. An important property of Eq. (13) is the existence of a manifold $U$ that orbits cannot cross. As shown in Appendix A, Eq. (13) is singular for the manifold determined by the condition

$$ \Gamma_R(\xi) = 0, \qquad (16) $$

with $\Gamma_R$ given by Eq. (A41). The role of this set is similar to the sonic circle found in the Hall-MHD model [24].

We first note that the upstream state $\xi_0 = [1 \ 0 \ 0 \ 0 \ 1]^T$ is an equilibrium state of Eq. (13) because it satisfies $\mathbf{f}(\xi_0) = 0$. Another interesting element is the stable (unstable) manifold $W^s(W^u)$ of $\xi_0$, which is the set of forward (backward) in $\hat{x}$ trajectories that terminate at $\xi_0$. Since solitary waves are localized structures that approach upstream and downstream to $\xi_0$, i.e., $\xi \to \xi_0$ as $\hat{x} \to \pm\infty$, these special solutions belong to the intersections of the stable and the unstable manifolds of $\xi_0$. They are called homoclinic orbits. As explained below, their organization in parameter space depends on (1) the dimension of the phase space, (2) the occurrence of symmetries, and (3) the dimensions of the stable and the unstable manifolds of $\xi_0$. These three topics are discussed below.

In the particular case of Eq. (13), the dimension of the phase space, given by the number of components of $\xi$, is five. However, as shown in Appendix A, there is a function $H(\xi)$ that satisfies $dH/d\hat{x} = 0$ [21,24], i.e., it is conserved. As a consequence, the effective dimension of our system is four. Regarding symmetries, one readily verifies that Eq. (13) is reversible because it admits the involution $G\mathbf{f}(\xi) = -\mathbf{f}(G\xi)$ with

$$ G: (u_x, u_y, u_z, b_y, b_z) \to (u_x, -u_y, u_z, -b_y, b_z). \qquad (17) $$

The subspace $S: u_y = b_y = 0$, a key element for the later computation of the solitary waves, is called the symmetric section of the reversibility. Interestingly, a symmetric solitary wave exists if the unstable manifold of $\xi_0$ intersects the symmetric section at a given point. The reason is that, by the reversibility property, $W^s$ should also intersect $S$ at the same point.

On the other hand, the tangent spaces of $W^s$ and $W^u$ have the same dimensions as the stable and unstable spaces of the linearization of $\mathbf{f}$ at $\xi_0$. Substituting $\xi = \xi_0 + \xi_1$ in Eq. (13) with $\xi_1$ a small perturbation and dropping higher order terms yield

$$ \frac{d\xi_1}{d\hat{x}} \approx \bar{\mathbf{J}}|_{\xi_0} \xi_1, \qquad (18) $$
<figure><img src="image_3.png" /><figcaption>FIG. 1. Characteristic velocities (left) and domains of stability of $\xi_0$ in the $C/V_A-\theta$ plane (right). The regions are (1) saddle-center, (2) focus-focus, (3) saddle-saddle, (4) and center-center. The parameters used for panels (a) and (b) [(c) and (d)] correspond to case 1 (2) in Table I. The green dotted lines in panels (a) and (c) are the velocities making $p_1 = 0$ in Eq. (19).</figcaption></figure>

where $\bar{J}|_{\xi_0}$ is the Jacobian matrix of $f$ at $\xi_0$. If we now assume $\xi_1(\hat{x}) = \hat{\xi}_1 e^{\lambda \hat{x}}$, the condition for nontrivial $\hat{\xi}_1$ is det($\bar{J}|_{\xi_0} - \lambda \bar{I}$) = 0. Such a condition gives $\lambda = 0$, which is a consequence of the invariant $H$, and the characteristic equation with a biquadratic structure that reflects the involution given by Eq. (17)

$$p_2\lambda^4 + p_1\lambda^2 + p_0 = 0, \quad (19)$$

where $p_2$, $p_1$, and $p_0$ are certain constants that just depend on the five parameters of Eq. (15) (one can find their explicit forms in Ref. [24]). These coefficients contain important information that will help us to connect the mathematical results with the physics of the solitary waves. Coefficient $p_0$ vanishes when the propagation velocity $C$ coincides with one of the nondispersive MHD velocities, i.e., the system obtained after neglecting the Hall and the FLR terms. The MHD velocities are the fast ($V_{fast}$) and the slow ($V_{slow}$) magnetosonic velocities and the firehose velocity ($V_F$), which reduces to the intermediate or shear Alfvén velocity in the case of isotropic pressure. Coefficients $p_1$ and $p_2$ vanish when the propagation velocity $C$ is equal to the acoustic velocity corrected with FLR effects ($V_s$) and the velocity $V_{FLR}$ defined by the condition $\Gamma_R(\xi_0) = 0$, respectively [24]. For the parameters of Table I, these velocities are plotted versus the propagation angle in Figs. 1(a) and 1(c).

The generic cases of the solutions of Eq. (19) are (1) saddle-center, $\lambda_{1,2} = \pm\kappa$ and $\lambda_{3,4} = \pm i\omega$, (2) saddle-saddle, $\lambda_{1,2} = \pm\kappa_1$ and $\lambda_{3,4} = \pm\kappa_2$, (3) focus-focus, $\lambda_{1,2} = \kappa \pm i\omega$ and $\lambda_{3,4} = -\kappa \pm i\omega$, and (4) center-center, $\lambda_{1,2} = \pm i\omega_1$ and $\lambda_{3,4} = \pm i\omega_2$. Figures 1(b) and 1(d) show the domains of stability of $\xi_0$ in the $C/V_A-\theta$ plane for cases (1) and (2) in Table I. Although this set of parameters yields to unstable solitary waves (see Sec. IV), they have been used throughout this work because they were used in Ref. [24]. Working with the same physical parameters eases the comparison of our

results and highlights the unique characteristics related with the existence of the solitary waves and their organization in parameter space. This particular case is also illustrative because, as shown in Fig. 1 the four stability regions of $\xi_0$, exist in the $C/V_A-\theta$ plane.

Taking into account that the effective dimension of the system is four and its reversible character, well-known theoretical results on the existence of homoclinic orbits can be directly applied to our case (one can find a review in Ref. [25]). To fix ideas, consider the situation with given $M_e, M_i$, and $a_p$ values and let us discuss the organization of solitary waves in the $M_A-\theta$ plane (as was already done in Ref. [24]). Unless very specific resonance conditions are fulfilled, no solitary wave occurs when $\xi_0$ is a center-center because such a point has no stable or unstable manifold and orbits cannot connect with it. For values of $M_A$ and $\theta$ making $\xi_0$ a saddle-center, the stable and unstable manifolds have dimension equal to one, and an homoclinic orbit exists if the two coincide, $W^s = W^u$. In general, the intersection of the one-dimensional manifold $W^u$ with the two-dimensional symmetric section is expected to occur for specific parameter values that form branches in the $M_A-\theta$ plane. For parameter values where $\xi_0$ is hyperbolic, i.e., saddle-saddle and focus-focus, $W^u$ has dimension two. The intersection of such a two-dimensional manifold with the two-dimensional symmetric section in a four-dimensional phase space is generic, and solitary waves are expected to exist in continuous regions in the $M_A-\theta$ plane. This is called a continuous spectrum. According to this discussion, we expect that in Figs. 1(b) and 1(d) we will find branches of solutions in region 1, a continuous spectrum in regions 2 and 3, and no solitary wave in region 4.

### III. FLR-HALL-MHD SOLITARY WAVES

#### A. Saddle-center domain

According to previous geometrical arguments, solitary waves are organized in branches within the saddle-center domain. These branches can be computed, and their existence proved rigorously, by using the following bisection algorithm (see details in Ref. [29]). For a given set of parameters, Eq. (13) is integrated with initial condition equal to

$$\xi(\hat{x} = 0) = \xi_0 + \epsilon\xi^u, \quad (20)$$

where $\epsilon$ is a small parameter ($10^{-9}$ in our calculations) and $\xi^u$ is the unstable eigenvector of $\bar{J}|_{\xi_0}$, i.e., the eigenvector with the positive and real eigenvalue. Such initial condition guarantees that the orbits leaves the equilibrium state from the linearization of the unstable manifold. The orbit is computed up to the value $\hat{x}^*$ that satisfies the condition $b_y(\hat{x}^*) = 0$. We then recorded the value of $u_y(\hat{x}^*)$. Such a procedure is repeated by covering a range of $C/V_A$ (or $\theta$) values, and we then construct a diagram with $u_y(\hat{x}^*)$ versus $C/V_A$ (or $\theta$). Each time a change of sign in $u_y(\hat{x}^*)$ occurs, it means that $u_y(\hat{x}^*)$ passes through zero and there is an orbit leaving the unstable manifold of $\xi_0$ and hitting the symmetric section. Therefore, a solitary wave exists.

Figure 2(a) shows $u_y(\hat{x}^*)$ versus $\theta$ for $C/V_A = 1$ and the parameters of case 1. The zeros of $u_y(\hat{x}^*)$ have been highlighted by plotting the absolute value of $u_y(\hat{x}^*)$ in logarithmic scale and using blue crosses and red dots for positive and negative
<figure><img src="image_4.png" /><figcaption>FIG. 2. Panels (a) and (b) show $u_y(\hat{x}^*)$ versus $\theta$ for $C/v_A = 1$, and $u_y(\hat{x}^*)$ versus $C/v_A$ for $\theta = 75^\circ$, respectively. Other parameters correspond to case 1 in Table I. Positive (negative) values of $u_y(\hat{x}^*)$ are denoted with blue crosses (red dots).</figcaption></figure>

values of $u_y(\hat{x}^*)$, respectively. Clearly, solitary waves exist for $\theta \approx 88.15^\circ$ and $\theta = 85.45^\circ$ (see inset) and $\theta = 48.7^\circ$. For $\theta < 47^\circ$, where the crosses and dots are mixed and do not follow a smooth curve, we cannot guarantee (neither rule out) the existence of branches. The reason is that the values of $u_y(\hat{x}^*)$ are very small and they fall below our numerical error, which is a combination of factors including the finite value of $\epsilon$, the integration error, and the finite precision arithmetic (double precision floating-point format used here).

Figure 2(b) shows a similar diagram but varying $C/v_A$ for $\theta = 75^\circ$. For $C/v_A > 1.15$ we find again a parametric region where $u_y(\hat{x}^*)$ is smaller than our error. For lower propagation velocity, as $C/v_A$ decreases, one first finds a wave with $C/v_A \approx 0.73$, a gap, a velocity range with many waves, a second gap, and another region with several waves (see inset). The gaps appear because, when launching an integration along the unstable manifold, the orbit hit the singular manifold $\Gamma_R(\xi) = 0$. An intensive parametric survey constructing diagrams such as the ones in Fig. 2 allowed us to present the branches of solitary waves in the two saddle-center regions of the $C/v_A-\theta$ plane (see Fig. 3). In the saddle-center region delimited by the firehose and the fast magnetosonic velocities there are many branches of solutions, especially close to $V_F$. For large propagation angles, it is even possible to find solitary waves with propagation velocities larger than the Alfvén velocity.

The limiting factor in overcoming the numerical problems in the calculations presented here is the use of finite precision arithmetic. This is indicated by the fact that further reducing $\epsilon$ or increasing the accuracy of the integration does not lead to resolution of solitary wave branches for $\theta < 50^\circ$ in Fig. 3. In order to progress further we implemented our method in the computer algebra system Mathematica, taking advantage of its arbitrary precision capabilities. Using 30 digits of working precision and error tolerance of 20 digits in the integrator and taking $\epsilon = 10^{-15}$ allow us to resolve branches of solitary waves in this problematic regime. We show one example

<figure><img src="image_5.png" /><figcaption>FIG. 3. Branches of solitary waves in the $C/V_A-\theta$ plane for case 1.</figcaption></figure>

branch as a dashed line in Fig. 3. This branch is tracked until numerical errors once again prevent us from isolating solitary wave solutions. Even though a further increase in precision could help proceed towards lower values of $\theta$, the computations quickly become very expensive, and we do not pursue an exhaustive determination of branches. The main point we illustrate here is that the difficulties in locating solitary waves for smaller angles are indeed numerical and can be overcome by increasing the precision of the calculations.

We note that works on electromagnetic solitary waves in relativistic plasma (laser-plasma interaction framework) have encountered similar difficulties. In particular, waves in regions such as the one shown in Fig. 2(a) with $\theta < 47^\circ$, where the residual value of $u_y(\hat{x}^*)$ is very small, were erroneously taken as true waves with a continuous spectrum in early works. It was later shown that they should be organized in branches in the saddle-center domain and the claimed waves were numerical artifacts; see Ref. [29] and references therein. Similarly, the FLR-Hall-MHD solitary wave presented in Fig. 2 of Ref. [24] is not a true localized structure because (1) the values of $\theta$ and $C/V_A$ were selected without looking for a branch (a relation between $\theta$ and $C/V_A$) and (2) the author found a value of $u_y(\hat{x}^*) \approx 10^{-7}$ within this (numerically difficult) parametric domain.

Figure 3 also shows that many solitary waves can exist with propagation velocity covering a broad range between the sonic and the fast magnetosonic velocities (and not only close to $V_{fast}$ as concluded in Ref. [24]). In the Hall-MHD model, these waves are of a dark type and were termed the fast magnetosonic family [13]. As shown below, the solitary waves found in the FLR-Hall-MHD model are also dark for that regime. Above point (d) there is a blank region because the orbits started along the unstable manifold hit the surface $\Gamma_R = 0$.

In order to illustrate the different types of solitary waves, we selected six cases in Fig. 3 and labeled them with letters from (a) to (f). For all of them, we plotted the velocity components, the modulus of the magnetic field normal to the propagation direction $b$, and the magnetic hodograph (see
<figure><img src="image_6.png" /><figcaption>FIG. 4. Solitary waves named (a) and (b) in Fig. 3.</figcaption></figure>

Figs. 4–6). In the latter, we denoted with an arrow the increasing direction of $\hat{x}$ and kept the same scale for both axes to ease the interpretation of the wave polarization. Since all the selected waves have $u_x > 1$ at $\hat{x} = 0$, the relation $\rho/\rho_0 = 1/u_x$ indicates that the densities at the center of the structures are lower than the background value.

Solitary waves (a)–(e) belong to the saddle-center domain with propagation velocities between the firehose and the fast magnetosonic velocities. Waves (a) and (b) are dark solitary waves, i.e., the magnetic field exhibits a minimum at the center of the structure (see Fig. 4). Wave (a), which has a larger propagation angle and velocity, exhibits a much lower depression of the magnetic field, and its polarization is more linear as compared with wave (b). Waves (c), (d), and (e), shown in Figs 5 and 6, have been selected to illustrate the set of branches that populate the central region of Fig. 3. For a given propagation angle, for instance, $\theta = 70^{\circ}$ in cases (c) and (d), the solitary waves develop more and more oscillations as the velocity decreases. The polarization is almost circular. Wave (e), which propagates almost normal to the ambient

<figure><img src="image_7.png" /><figcaption>FIG. 6. Solitary waves named (e) and (f) in Fig. 3.</figcaption></figure>

magnetic field, has magnetic field variations of order unity, but much stronger changes on the normalized velocity components (up to 40 times the propagation velocity). The central core of the solitary wave is also complex and involves several peaks.

The saddle-center domain enclosed by the slow magnetosonic and the FLR velocities is particularly interesting from a physical point of view. No wave was found numerically in Ref. [24] for this domain with the FLR-Hall-MHD model. Moreover, no solitary wave exists in the Hall-MHD model because the upstream state is of the center type. As shown in Fig. 3, the FLR effect opens new possibilities because a branch of solutions occurs with propagation angle between $70^{\circ} \le \theta \le 74^{\circ}$. The solitary wave named (f), which is an example of such a branch, shows that for this domain the solitary waves are of the bright type, i.e., they exhibit a maximum of the magnetic field at its center. This slow family presents small (large) modulations of the magnetic (velocity) field components (note the different scale of the left and right axes in the $u_x$-$\hat{x}$ and $b$-$\hat{x}$ diagrams).

## B. Saddle-saddle and focus-focus domains

In the saddle-saddle and focus-focus cases, since the dimension of the unstable manifold is two and the solitary waves have a continuous spectrum, the algorithm should be modified slightly. For given parameter values, the initial condition in the saddle-saddle and focus-focus cases are

$$ \xi(\hat{x} = 0) = \xi_0 + \epsilon(\cos \varphi \xi_1^u + \sin \varphi \xi_2^u), \quad (21) $$

$$ \xi(\hat{x} = 0) = \xi_0 + \epsilon \mathrm{Re}(e^{i\varphi} \xi^u), \quad (22) $$

with $\xi_1^u$ and $\xi_2^u$ the eigenvectors with positive eigenvalues for the saddle-saddle case, and $\xi^u$ any of the two eigenvectors with eigenvalues having a positive real part for the focus-focus case. Angle $\varphi$ is a numerical parameter that controls the position of the initial condition in the linearization of the unstable manifold. The numerical scheme is similar to the saddle-center case, but now we need to look for the change of sign of $u_y(\hat{x}^*)$ as a function of $\varphi$.

<figure><img src="image_8.png" /><figcaption>FIG. 5. Solitary waves named (c) and (d) in Fig. 3.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. | $u_y(\hat{x}^*)|$ versus $\varphi$ diagram (a)] and example of solitary wave in the focus-focus domain (b)–(d).</figcaption></figure>

In order to illustrate this case, we now present some results for the parameter values of case (2) in Table I, $\theta = 70^{\circ}$ and $C/V_A = 0.52$. We set the numerical parameter $\epsilon = 10^{-9}$ and computed the orbits of Eq. (13) with initial conditions given by Eq. (22) and $\varphi$ from $96^{\circ}$ to $112^{\circ}$. For each of them, the value of $u_y$ at the intersection with the symmetric section, $u_y(\hat{x}^*)$, was computed and presented in a $|u_y(\hat{x}^*)|$ versus $\varphi$ diagram [see Fig. 7(a)]. Similarly to the previous section, blue crosses and red dots were used to denote positive and negative values of $u_y(\hat{x}^*)$ and highlight the changes of signs and locate the existence of a solitary waves. For instance, a solitary wave exists for $\varphi \approx 100.3^{\circ}$, and its structure is given in Figs. 7(b)–7(d). Interestingly, it can be proved that the existence of one solitary wave for a given value of the physical parameters implies the existence of infinitely many others if the system is reversible and the upstream state is a focus-focus [30]. Such a theoretical result, which was demonstrated earlier for conservative systems [31], is a consequence of the spiralling linear dynamics due to the complex eigenvalue and the additional orbits are like copies of the original one but with extra oscillations. Several of these extra orbits can be identified in Fig. 7(a).

## IV. STABILITY OF SOLITARY WAVES

The double adiabatic model for the pressures has been very useful in Secs. II and III because the dimension of the phase space is low and well-known techniques from the theory of the dynamical systems can be used. The organization of the solitary waves in the propagation angle-velocity plane was studied and their main physical features, such as polarizations and structure, found. However, just the existence of these solutions in the FLR-Hall-MHD model does not guarantee their physical relevance. The observation of these localized structures in real scenarios, such as the solar wind, is also linked to the concept of stability. For this reason, we now present some numerical simulations of the nonstationary FLR-Hall-MHD [Eqs. (B1)–(B3)] equations initialized with the solitary waves found in Sec. III. As we will see, the propagation of certain solutions is limited by an instability

<figure><img src="image_10.png" /><figcaption>FIG. 8. Panels (a) and (b) show the maximum growth rate $\gamma$ in the $k-\theta$ plane for cases 1 and 2 in Table I, respectively.</figcaption></figure>

of the background plasma state if the equations are closed with the double adiabatic pressure model (see Sec. IV A). However, we show that such solutions are useful because they are robust and can propagate indefinitely if a more accurate framework, i.e., with evolution equations for the pressures, is implemented.

### A. Dispersion relation of the FLR-Hall MHD system

As a preliminary step, we first analyzed the linear stability, i.e., dispersion relation, of the background plasma state with a double adiabatic pressure model. Since the solitary waves satisfy $\xi \to \xi_0$ as $\hat{x} \to \pm\infty$, a necessary condition for their stability is the linear stability of the background plasma state. We analyze it by writing the fluid variables as

$$ \hat{\rho} = 1 + \hat{\rho}_1 e^{i(k\hat{x} - \omega\tau)}, \quad (23) $$

$$ \mathbf{u} = \mathbf{e}_x + (\hat{u}_{x1}\mathbf{e}_x + \hat{u}_{y1}\mathbf{e}_y + \hat{u}_{z1}\mathbf{e}_z)e^{i(k\hat{x} - \omega\tau)}, \quad (24) $$

$$ \hat{\mathbf{B}} = \frac{\mathbf{e}_x}{\tan \theta} + \mathbf{e}_z + (\hat{b}_{y1}\mathbf{e}_y + \hat{b}_{z1}\mathbf{e}_z)e^{i(k\hat{x} - \omega\tau)}, \quad (25) $$

where $\hat{\rho} = \rho/\rho_0$, $\hat{B} = B/B_0 \sin \theta$, and $k$ and $\omega$ represent the normalized wave vector and frequency, respectively, of the small perturbations denoted with subscript 1. Substituting these expansions in Eqs. (B1)–(B3) and retaining only first-order terms yields the homogeneous linear system

$$ \tilde{D}(\omega, k)\hat{\eta} = 0, \quad (26) $$

with $\hat{\eta} = [\hat{\rho}_1 \hat{u}_{x1} \hat{u}_{y1} \hat{u}_{z1} \hat{b}_{y1} \hat{b}_{z1}]$. For given values of $M_e, M_i, \theta$, and $a_p$, the compatibility condition $\det(\tilde{D}) = 0$ gives the dispersion relation $\omega = \omega(k)$. The background state is unstable if $\omega$ is imaginary, and we can compute the growth rate as $\gamma = \max[\tilde{\zeta}(\omega)]$.

Figures 8(a) and 8(b) show the value of the growth rate in the $k-\theta$ plane for cases 1 and 2 in Table I. For this particular set of parameters we conclude that only low-angle propagation waves could be stable. By comparing this diagram with Fig. 3, one finds that the tails of all the waves computed at $\theta > 50^{\circ}$
<figure><img src="image_11.png" /><figcaption>FIG. 9. Evolutions of some solitary waves belonging to the saddle-center regime. Panel (a) shows an example of unstable core for $\theta = 80^{\circ}$ and $C/v_A = 0.745715$. Panel (b) shows an example of unstable background for $\theta = 30.415^{\circ}$ and $C/v_A = 0.9$. Panel (c) shows an example of robust solution, using Eqs. (27) and (28) for the pressures with the same parameters as panel (b).</figcaption></figure>

are unstable. The use of extended precision allowed us to locate solitary wave solutions at lower propagation angles, for which the maximum growth rate of the perturbations to the background state tends to zero.

## B. Numerical simulations

This section studies the stability of two solitary waves with parameters given by case 1 in Table I. In both cases, the waves belong to the parameter regime where $\xi_0$ is a saddle-center and a relation between $\theta$ and $C/V_A$ holds (branches of solutions). They were used as initial conditions in Eqs. (B1)–(B3), and their evolutions were found by integrating the equations numerically with a spectral method, following Ref. [32] (some details on the numerical method are found in Appendix B).

The evolution of the spatial profile of $u_x$ for the first wave, which has $\theta = 80^{\circ}$ and $C/V_A = 0.745715$, is shown in Fig. 9(a). According to Fig. 8 the background plasma state is unstable for such a high propagation angle. However, as shown in Fig. 9, the core of the solitary wave is unstable, and the solitary wave is destroyed even before the instability at the tail would be developed. The behavior of the second wave, having velocity $C/V_A = 0.9$ and a propagation angle $\theta = 30.415^{\circ}$, is totally different [see Fig. 9(b)]. For this case, the core of the solitary wave is stable, and the instability happens at the tail. The results of the simulation, i.e., the most

unstable wave vector and the growth rate, are consistent with the analysis of Sec. IV A. Interestingly, although the wave is unstable, the core of this solitary wave is quite robust and survives a time longer than $250 M_A \cos \theta / \Omega_{ci,0} \approx 266 \Omega_{ci,0}^{-1}$.

The simulation results for the second wave show that the instability may come from the unstable character of the background plasma state in the framework of the FLR-Hall-MHD system closed with a double adiabatic pressure model. For this reason, we investigated a bit further the stability of the second wave by repeating the simulations but now using the dynamic equations for the pressures. Following Ref. [26] (a short discussion is found in Appendix C), we write

$$ \frac{\partial P_{\parallel}}{\partial \tau} + \hat{\nabla} \cdot (P_{\parallel} \mathbf{u}) + 2P_{\parallel} \mathbf{e}_b \cdot \hat{\nabla} \mathbf{u} \cdot \mathbf{e}_b + \frac{1}{a_p} \left[ (\bar{\mathbf{\Pi}} \cdot \hat{\nabla} \mathbf{u})^S : \bar{\tau} - \bar{\mathbf{\Pi}} : \frac{d\bar{\tau}}{d\tau} \right] = 0, \quad (27) $$

$$ \frac{\partial P_{\perp}}{\partial \tau} + \hat{\nabla} \cdot (P_{\perp} \mathbf{u}) + P_{\perp} \hat{\nabla} \cdot \mathbf{u} - P_{\perp} \mathbf{e}_b \cdot \hat{\nabla} \mathbf{u} \cdot \mathbf{e}_b + \frac{1}{2} \left[ (\bar{\mathbf{\Pi}} \cdot \hat{\nabla} \mathbf{u})^S : \bar{\mathbf{I}} - (\bar{\mathbf{\Pi}} \cdot \hat{\nabla} \mathbf{u})^S : \bar{\tau} + \bar{\mathbf{\Pi}} : \frac{d\bar{\tau}}{d\tau} \right] = 0, \quad (28) $$

where $\hat{\nabla} = \partial/\partial \hat{x}$, $\bar{\tau} = \mathbf{e}_b \cdot \mathbf{e}_b$ and $\bar{\mathbf{\Pi}}$ is the nondimensional ion FLR pressure tensor, given by Eqs. (6)–(8) scaled with $p_{\perp 0}$. The superscript S means that the tensor between parentheses is symmetrized by the addition of its transpose. Unlike the double adiabatic model, Eqs. (B1)–(B3) and Eqs. (27) and (28) conserve the energy

$$ E = \int_{-\infty}^{\infty} \left[ \frac{1}{2} \hat{\rho} |\mathbf{u}|^2 + \frac{1}{2} M_A \sin^2 \theta b^2 + M_e \hat{\rho} \log \hat{\rho} + M_i \left( P_{\perp} + \frac{1}{2} a_p P_{\parallel} \right) \right] d\hat{x}. \quad (29) $$

As shown in Fig. 9(c), the substitution of the crude double adiabatic approximation by the pressure evolution equation suppresses the instability of the wave. Since the initial condition is not an exact solution of the FLR-Hall-MHD model closed with Eqs. (27) and (28), the wave is distorted slightly, but it still propagates for times longer than $500 M_A \cos \theta / \Omega_{ci,0}$ while keeping its original shape. The time integration was stopped at $500 M_A \cos \theta / \Omega_{ci,0}$, but the simulation was still stable. The shape of the wave at this time is practically identical to the given initial condition, and it propagates with the speed of the used reference frame (the wave does not drift). These results suggest that this particular solitary wave computed with the double adiabatic pressure model is very close to an exact solution of the system with dynamical pressure equations. Similar simulations performed at $\theta = 50^{\circ}$ showed less robust behavior, greater level of deformation, and drift leftwards of the domain. These interesting results justify the approach followed in this work, which provided useful initial conditions for solitary waves obtained from the FLR-Hall-MHD closed with Eqs. (27) and (28).

## V. CONCLUSIONS

The existence of low-frequency solitary waves in magnetized plasmas is firmly supported by space observations. For this reason, the knowledge of the physical properties of these
structures, including amplitudes, spatial profiles of the fluid and electromagnetic fields, and polarizations, are relevant. The existence of possible relations linking physical parameters, such as propagation angle of the wave with respect to the ambient magnetic field $\theta$ and the propagation velocity $C$, are also important because they can be helpful in interpreting the experimental data. The answers to most of these interesting questions can be obtained by analyzing the dynamical system obtained from the double adiabatic FLR-Hall-MHD model after assuming the one-dimensional traveling wave ansatz.

First, solitary waves can exist if the background plasma state, which appears in the dynamical system as an equilibrium state, is not a center-center. Moreover, using simple geometrical arguments based on the effective dimension of the dynamical system and its reversible character, the organization of the waves in the $\theta-C$ plane can be anticipated even before computing them. If the background plasma state is a saddle-saddle or a focus-focus the spectrum of the waves is continuous, and, in the case where it is a saddle-center, they are organized in branches [relations of the type $C = C(\theta)$]. The numerical scheme (bisection method) presented in this work can be used to find solitary waves in any of these regions and proves their existence rigorously. Abundant solitary waves, including dark and bright waves and some of them belonging to regions where to the best of our knowledge they were not found before, were computed. The structures of seven waves were presented in detail, and some differences with respect to the Hall-MHD case (without FLR effects) were highlighted. Nevertheless, deeper parametric analysis is necessary to construct a more complete picture about the properties and organization of the waves in parameter space. For instance, the fact that we did not find waves with banana-like polarization, a very peculiar signature observed in the solar wind and in more simple theoretical models, does not preclude their existence in this FLR-Hall-MHD model. Another topic that could be investigated in future works is the analysis of the existence of the so-called quasisolitons, i.e., a more general class of solutions that would contain the branches of solutions found in this work as a particular subclass.

Regarding the stability of the waves, a linear analysis within the framework of the double adiabatic FLR-Hall-MHD model shows that the background plasma state is unstable for the parameters under consideration. Since the tails of the waves approach such a state at plus and minus infinity, they are also unstable. However, some numerical simulations indicate that the core can be stable. It is remarkable that closing the FLR-Hall-MHD system by the double adiabatic pressure model is not only problematic because of the instability of the background plasma state (even if isotropic with $a_p = 1$) but also because energy conservation is not guaranteed. Even without FLR effects, i.e., for the Hall-MHD model, an energy conservation theorem is not known. From this perspective, it is not surprising that these nonphysical features disappear if the double adiabatic model is substituted by evolution equations for the pressures, which conserve the energy and allow stable states of the background plasma. For that case, it has been shown that some of the solitary waves are robust. This result opens the interesting problem about the computation of exact solitary waves in the framework of the FLR-Hall-MHD with evolution equations for the pressures. Such a study, which is

beyond the scope of the present work, is challenging because the dimension of the phase space of the dynamical system would be larger and several arguments used in this work should be revised.

## ACKNOWLEDGMENTS

G.S.A. is supported by the Ministerio de Economía y Competitividad of Spain under Grant No. RYC-2014-15357. E.S. was supported by the Swedish Research Council, Grant No. 2016-05012, and by the Knut and Alice Wallenberg Foundation.

## APPENDIX A: DYNAMICAL SYSTEM

This Appendix follows Ref. [24] to find the explicit form of vector $\mathbf{f}$ in Eq. (13). For convenience, we split such a column vector as $\mathbf{f} = [\mathbf{f}_u \ \mathbf{f}_b]$, with $\mathbf{f}_u$ and $\mathbf{f}_b$ the three- and two-dimensional vector flows governing the dynamics of $\mathbf{u}$ and $\mathbf{b}$, respectively. Equations (1)–(3) are particularized to one-dimensional ($\partial/\partial y = \partial/\partial y = 0$) and steady ($\partial/\partial t = 0$) solutions. Equation (1) then becomes

$$ \frac{\rho}{\rho_0} = \frac{1}{u_x}. \qquad (A1) $$

After defining the new variables $P_{\parallel}(u, b^2) \equiv p_{\perp}/p_{\perp 0}$ and $P_{\parallel}(u, b^2) \equiv p_{\parallel}/p_{\parallel 0}$, the equations of state (9) and (10) read

$$ P_{\parallel}(u_x, b^2) = \frac{1}{\hat{b}^2 u_x^3}, \quad P_{\perp}(u_x, b^2) = \frac{\hat{b}}{u_x}, \qquad (A2) $$

where we introduced the dimensionless quantities $b^2 = b_y^2 + b_z^2$ and $\hat{b}^2 = (B/B_0)^2 = \cos^2 \theta + b^2 \sin^2 \theta$.

The component of Eq. (3) along the propagation direction $x$ gives $B_x = B_0 \cos \theta$. In the transverse direction one finds

$$ \frac{d}{d\hat{x}} \left[ u_x b_z \sin \theta - u_z \cos \theta + u_x \sin \theta \frac{db_y}{d\hat{x}} \right] = 0, \qquad (A3) $$

$$ \frac{d}{d\hat{x}} \left[ u_x b_y \sin \theta - u_y \cos \theta - u_x \sin \theta \frac{db_z}{d\hat{x}} \right] = 0. \qquad (A4) $$

Using the plasma conditions upstream, this set of equations are integrated to find the two-dimensional flow

$$ f_b = \left( \begin{array}{c} \frac{u_x}{u_x} \frac{\cos \theta}{\sin \theta} - b_z + \frac{1}{u_x} \\ - \frac{u_y}{u_x} \frac{\cos \theta}{\sin \theta} + b_y \end{array} \right). \qquad (A5) $$

This flow coincides with Eq. (17) in Ref. [24], except for the term $1/u_x$ in the first row of Eq. (A5).

Following a similar procedure, Eq. (2) gives

$$ \bar{\mathbf{A}} \cdot \frac{d\mathbf{u}}{d\hat{x}} + \mathbf{F}(\mathbf{u}, \mathbf{b}) = 0, \qquad (A6) $$

where we introduced the flow $\mathbf{F} = (F_x \mathbf{e}_x + F_y \mathbf{e}_y + F_z \mathbf{e}_z)/\delta$ with

$$ F_x = u_x - 1 + P(u_x, b^2) + \frac{1}{2} M_A \sin^2 \theta (b^2 - 1), \qquad (A7) $$

$$ F_y = u_y + \chi(u_x, b^2) \cos \theta \sin \theta b_y, \qquad (A8) $$
$$F_z = u_z + [\chi(u_x, b^2)b_z - \chi(1, 1)] \cos \theta \sin \theta, \quad (A9)$$

and the auxiliary functions

$$\delta = \frac{M_i P_{\perp}(u_x, b^2)}{M_A \hat{b} \cos \theta}, \quad (A10)$$

$$P(u_x, b^2) = M_e \left(\frac{1}{u_x} - 1\right) + M_i \left\{P_{\perp}(u_x, b^2) - 1 + [a_p P_{\parallel}(u_x, b^2) - P_{\perp}(u_x, b^2)] \frac{\cos^2 \theta}{\hat{b}^2} - (a_p - 1) \cos^2 \theta \right\}, \quad (A12)$$

$$\chi(u_x, b^2) = M_i [a_p P_{\parallel}(u, b^2) - P_{\perp}(u, b^2)] \frac{1}{\hat{b}^2} - M_A. \quad (A13)$$

Factor $1/\hat{b}$ appearing in Eq. (A10), which comes from the fact that $\Omega_{ci}$ in Eqs. (6)–(8) is the local ion gyro frequency, is missed in Ref. [24]. Tensor $\bar{\mathbf{A}}$ in Equation (A6) is

$$\bar{\mathbf{A}} = \bar{\mathbf{I}} \times \mathbf{r} - 2\epsilon \hat{\mathbf{b}}_{\parallel} \times \hat{\mathbf{b}}_{\perp}, \quad (A14)$$

with $\hat{\mathbf{b}}_{\parallel} = \hat{b}_x \mathbf{e}_x$, $\hat{\mathbf{b}}_{\perp} = \hat{b}_y \mathbf{e}_y + \hat{b}_z \mathbf{e}_z$, $\hat{\mathbf{b}} = \hat{\mathbf{b}}_{\parallel} + \hat{\mathbf{b}}_{\perp}$, $\hat{b}_x = B_x/B = \cos \theta/\hat{b}$, and $\hat{b}_{y,z} = B_{y,z}/B = b_{y,z} \sin \theta/\hat{b}$, and

$$\mathbf{r} = -r_{\parallel} \hat{\mathbf{b}}_{\parallel} + r_{\perp} \hat{\mathbf{b}}_{\perp}, \quad (A15)$$

$$r_{\parallel} = \frac{1}{2}(1 - 3\hat{b}_{\parallel}^2) + 2\epsilon\hat{b}_{\parallel}^2, \quad (A16)$$

$$r_{\perp} = \frac{1}{2}(1 + 3\hat{b}_{\parallel}^2) - 2\epsilon\hat{b}_{\parallel}^2, \quad (A17)$$

$$\epsilon = (p_{\perp} - p_{\parallel})/p_{\perp} = 1 - a_p P_{\parallel}/P_{\perp}. \quad (A18)$$

Following Ref. [24], tensor $\bar{\mathbf{A}}$ will be referred as the 1D FLR tensor. We mention that a plus sign (instead a minus) was written in the second term of Eq. (A14) in Ref. [24].

## 1. Singularity of the tensor $\bar{\mathbf{A}}$ and invariant manifold

As pointed out in Refs. [21] and [24], tensor $\bar{\mathbf{A}}$ is singular:

$$\mathbf{L} \cdot \bar{\mathbf{A}} = 0, \quad (A19)$$

$$\bar{\mathbf{A}} \cdot \mathbf{R} = 0, \quad (A20)$$

and left and right null vectors are given by

$$\mathbf{L} = \frac{1}{\mu}(\mathbf{r} + 2\epsilon\hat{b}_{\perp}^2\hat{\mathbf{b}}_{\parallel} - 2\epsilon\hat{b}_{\parallel}^2\hat{\mathbf{b}}_{\perp}), \quad (A21)$$

$$\mathbf{R} = \frac{1}{\mu}\mathbf{r}. \quad (A22)$$

After imposing the condition $\mathbf{L} \cdot \mathbf{R} = 1$, the arbitrary constant $\mu$ is

$$\mu^2 = r_{\perp}^2 \hat{b}_{\perp}^2 + r_{\parallel}^2 \hat{b}_{\parallel}^2 - \gamma, \quad (A23)$$

with $\gamma = 2\epsilon\hat{b}_{\perp}^2\hat{b}_{\parallel}^2$. A direct result of the singular character of $\bar{\mathbf{A}}$ is the constraint

$$H(\xi) = \mathbf{L} \cdot \mathbf{F} = 0, \quad (A24)$$

which is easily obtained by dotting Eq. (A6) from the left with $\mathbf{L}$. Therefore, any orbit of the five-dimensional phase

space of $\xi$ in Eq. (13) should lie in the four-dimensional surface defined by the constraint $H$. As a consequence, the effective dimension of the system is four. Although a detailed derivation on how Eq. (13) can be obtained from (A6) was given in Ref. [24], we summarize below the most important calculations because some small discrepancies were found.

Besides the zero eigenvalue, the 1D FLR tensor $\bar{\mathbf{A}}$ has imaginary eigenvalues $\pm i\mu$, where $\mu$ is given by Eq. (A23). It can be shown [24] that vectors

$$\mathbf{S} = \frac{1}{\mu \hat{b}_{\perp} \hat{b}_{\parallel}} [(\hat{b}_{\perp}^2 r_{\perp} - \gamma) \hat{\mathbf{b}}_{\parallel} + (\hat{b}_{\parallel}^2 r_{\parallel} - \gamma) \hat{\mathbf{b}}_{\perp}], \quad (A25)$$

$$\mathbf{T} = \frac{1}{\hat{b}_{\perp} \hat{b}_{\parallel}} (\hat{\mathbf{b}}_{\parallel} \times \hat{\mathbf{b}}_{\perp}), \quad (A26)$$

$$\mathbf{M} = \frac{1}{\mu \hat{b}_{\perp} \hat{b}_{\parallel}} (\hat{b}_{\perp}^2 r_{\perp} \hat{\mathbf{b}}_{\parallel} + \hat{b}_{\parallel}^2 r_{\parallel} \hat{\mathbf{b}}_{\perp}), \quad (A27)$$

$$\mathbf{N} = \mathbf{T} \quad (A28)$$

satisfy the relations $\bar{\mathbf{A}} \cdot \mathbf{S} = -\mu \mathbf{T}$, $\bar{\mathbf{A}} \cdot \mathbf{T} = \mu \mathbf{S}$, $\mathbf{M} \cdot \bar{\mathbf{A}} = \mu \mathbf{N}$, and $\mathbf{N} \cdot \bar{\mathbf{A}} = -\mu \mathbf{M}$. One also readily finds that the following orthogonality and normalization conditions hold:

$$\mathbf{L} \cdot \mathbf{R} = \mathbf{M} \cdot \mathbf{S} = \mathbf{T} \cdot \mathbf{N} = 1, \quad (A29)$$

$$\mathbf{M} \cdot \mathbf{R} = \mathbf{N} \cdot \mathbf{R} = \mathbf{L} \cdot \mathbf{S} = \mathbf{L} \cdot \mathbf{T} = \mathbf{N} \cdot \mathbf{S} = \mathbf{M} \cdot \mathbf{T} = 0. \quad (A30)$$

The new base $\mathbf{R}$, $\mathbf{S}$, $\mathbf{T}$ will allow us to find the flow $\mathbf{f}_u$ in Eq. (13) from Eq. (A6). We first decompose $\mathbf{f}_u$ and $\mathbf{F}$ on that base and write

$$\mathbf{F} = F_R \mathbf{R} + F_S \mathbf{S} + F_T \mathbf{T}, \quad (A31)$$

$$\mathbf{f}_u = f_u \mathbf{R} + f_u \mathbf{S} + f_u \mathbf{T}. \quad (A32)$$

The dot product of Eq. (A31) by $\mathbf{L}$, $\mathbf{M}$, and $\mathbf{N}$ gives

$$F_R = \mathbf{L} \cdot \mathbf{F} = 0, \quad (A33)$$

$$F_S = \mathbf{M} \cdot \mathbf{F}, \quad (A34)$$

$$F_T = \mathbf{N} \cdot \mathbf{F}, \quad (A35)$$

where we used Eq. (A24) and the orthogonality and normalization conditions. We now find the components of $\mathbf{f}_u$ by first noting that Eqs. (13) and (A6) give

$$\bar{\mathbf{A}} \cdot \mathbf{f}_u = -\mathbf{F}. \quad (A36)$$

The substitution in Eq. (A36) of Eq. (A32) and the use of Eqs. (A25) and (A26) yield

$$f_u S = F_T / \mu, \quad (A37)$$

$$f_u T = -F_S / \mu. \quad (A38)$$

Finally, the component $f_u R$ is found from the constraint (A24). From such invariant, one finds

$$\frac{dH}{d\hat{x}} = \frac{\partial H}{\partial \mathbf{u}} \cdot \frac{d\mathbf{u}}{d\hat{x}} + \frac{\partial H}{\partial \mathbf{b}} \cdot \frac{d\mathbf{b}}{d\hat{x}} = \frac{\partial H}{\partial \mathbf{u}} \cdot \mathbf{f}_u + \frac{\partial H}{\partial \mathbf{b}} \cdot \mathbf{f}_b = 0, \quad (A39)$$
where we used Eq. (13). The component $f_{uR}$ then reads

$$f_{uR} = - \frac{\Gamma_S f_{uS} + \Gamma_b}{\Gamma_R}, \qquad (A40)$$

with

$$\Gamma_R = \frac{\partial H}{\partial \mathbf{u}} \cdot \mathbf{R}, \quad \Gamma_S = \frac{\partial H}{\partial \mathbf{u}} \cdot \mathbf{S}, \quad \Gamma_b = \frac{\partial H}{\partial \mathbf{b}} \cdot \mathbf{f}_b, \quad (A41)$$

and where we used that $\partial H/\partial \mathbf{u} \cdot \mathbf{T} = 0$ because (1) $\mathbf{L} \cdot \mathbf{T} = 0$, (2) the derivatives of $\mathbf{L}$ with respect to $u_y$ and $u_z$ are all zero and $\partial \mathbf{L}/\partial \mathbf{u} \cdot \mathbf{F}$ is along $\mathbf{e}_x$, and (3) as shown by Eq. (A26), $\mathbf{T}$ is perpendicular to $\mathbf{e}_x$. The analytical derivatives of $\partial H/\partial \mathbf{u}$ and $\partial H/\partial \mathbf{b}$ have been implemented in our code. Equations (A39) and (A40) have a different sign compared with the corresponding equations in Ref. [24].

The initial conditions used in this work are consistent with the constraint (A24) because $H(\xi_0) = 0$. The flow $f$ in Eq. (13) guarantees that the orbit $\xi(\hat{x})$ will lie in the manifold $H = 0$. As pointed out in Ref. [24], orbits cannot cross the set $U$ defined by $\Gamma_R(U) = 0$, which plays a similar role to the sonic circle in the Hall-MHD theory [13].

## APPENDIX B: ONE-DIMENSIONAL FLR-HALL MHD SPECTRAL CODE

For convenience, the simulations in Sec. IV used the same dimensionless variables as in the previous section, and $\hat{\rho} = \rho/\rho_0$ and the normalized time $\tau = v_{x0}t/\ell$. After substituting $\nabla = \partial/\partial x \mathbf{e}_x$, Eqs. (1)–(3) become

$$\frac{\partial \hat{\rho}}{\partial \tau} + \frac{\partial}{\partial \hat{x}}(\hat{\rho} u_x) = 0, \qquad (B1)$$

$$\frac{\partial}{\partial \tau}(\hat{\rho}\mathbf{u}) + \mathbf{e}_x \cdot \frac{\partial}{\partial \hat{x}} \left[ \hat{\rho}\mathbf{u}\mathbf{u} + M_e \hat{\rho} \bar{\mathbf{I}} + M_i \hat{\rho} \hat{\mathbf{P}}_i^{(0)} \right. \\ \left. + M_A \sin^2 \theta \left( \frac{1}{2} \hat{B}^2 \bar{\mathbf{I}} - \hat{B}\hat{\mathbf{B}} \right) + M_i \bar{\mathbf{\Pi}} \right] = 0, \qquad (B2)$$

$$\frac{\partial \hat{\mathbf{B}}}{\partial \tau} = \mathbf{e}_x \times \left[ \frac{\partial}{\partial \hat{x}} \left( \mathbf{u} \times \hat{\mathbf{B}} - \frac{1}{\hat{\rho}} \frac{\partial \mathbf{b}}{\partial \hat{x}} \right) \right], \qquad (B3)$$

with

$$\hat{\mathbf{P}}_i^{(0)} = \frac{a_p \hat{\rho}^3}{\hat{b}^2} \mathbf{e}_b \mathbf{e}_b + \hat{\rho} \hat{b} (\bar{\mathbf{I}} - \mathbf{e}_b \mathbf{e}_b), \qquad (B4)$$

$$\hat{\mathbf{B}} = \frac{\mathbf{e}_x}{\tan \theta} + \mathbf{b}. \qquad (B5)$$

Tensor $\bar{\mathbf{\Pi}}$ accounts for the FLR effect, and only its first row

$$\mathbf{e}_x \cdot \bar{\mathbf{\Pi}} = \frac{1}{M_A \cos \theta} \frac{1}{\hat{b}} \bar{\mathbf{M}} \cdot \frac{\partial \mathbf{u}}{\partial \hat{x}} \qquad (B6)$$

is needed, where $\bar{\mathbf{M}} = P_\perp \bar{\mathbf{A}}$. After using the double adiabatic equations, tensor $\bar{\mathbf{A}}$ is given by Eq. (A14) but with $\varepsilon$ now taking the form $\varepsilon = 1 - a_p \hat{\rho}^2/\hat{b}^3$. The components of $\bar{\mathbf{M}}$ read

$$M_{11} = 0, \qquad (B7)$$

$$M_{12} = -\frac{1}{2} \hat{b}_z [P_\perp + (8a_p P_\parallel - 5P_\perp) \hat{b}_x^2], \qquad (B8)$$

$$M_{13} = \frac{1}{2} \hat{b}_y [P_\perp + (8a_p P_\parallel - 5P_\perp) \hat{b}_x^2], \qquad (B9)$$

$$M_{21} = \frac{1}{2} \hat{b}_z [P_\perp + (4a_p P_\parallel - P_\perp) \hat{b}_x^2], \qquad (B10)$$

$$M_{22} = 2\hat{b}_x\hat{b}_y\hat{b}_z(P_\perp - a_p P_\parallel), \qquad (B11)$$

$$M_{23} = \frac{1}{4} \hat{b}_x [P_\perp (3 + \hat{b}_x^2 - 9\hat{b}_y^2 - \hat{b}_z^2) + 8a_p P_\parallel (\hat{b}_y^2 - \hat{b}_x^2)], \qquad (B12)$$

$$M_{31} = \frac{1}{2} \hat{b}_y [P_\perp (\hat{b}_x^2 - 1) - 4a_p P_\parallel \hat{b}_x^2], \qquad (B13)$$

$$M_{32} = -\frac{1}{4} \hat{b}_x [P_\perp (3 + \hat{b}_x^2 - \hat{b}_y^2 - 9\hat{b}_z^2) + 8a_p P_\parallel (\hat{b}_z^2 - \hat{b}_x^2)], \qquad (B14)$$

$$M_{33} = -2\hat{b}_x\hat{b}_y\hat{b}_z(P_\perp - a_p P_\parallel). \qquad (B15)$$

Equations (B7)–(B15) have been written in terms of $P_\perp$ and $P_\parallel$. This is convenient since several closures for the ion pressure are being used at different stages of the work.

Equations (B1)–(B3) have been integrated numerically with the spectral method (see, e.g., Ref. [32]). The size of the simulation box and the number of points of the mesh (after desaliasing by a factor two) were equal to 94.328 and 2048, respectively. A spectral cutoff is imposed on the spectrum at half the spectral domain.

## APPENDIX C: FLR WORK IN DYNAMICAL PRESSURE EQUATIONS

This Appendix provides explicit equations for the terms appearing in the right-hand side of Eqs. (27) and (28). Particularizing for $\hat{\mathbf{V}} = \partial/\partial \hat{x}$ and after some development, the FLR work terms in these equations read

$$(\bar{\mathbf{\Pi}} \cdot \hat{\mathbf{V}} \mathbf{u})^S : \bar{\mathbf{\tau}} = 2(\mathbf{e}_b \cdot \bar{\mathbf{\Pi}} \cdot \mathbf{e}_x) \left( \mathbf{e}_b \cdot \frac{\partial \mathbf{u}}{\partial \hat{x}} \right), \qquad (C1)$$

$$(\bar{\mathbf{\Pi}} \cdot \hat{\mathbf{V}} \mathbf{u})^S : \bar{\mathbf{I}} = 2 \frac{\partial \mathbf{u}}{\partial \hat{x}} \cdot \bar{\mathbf{\Pi}} \cdot \mathbf{e}_x, \qquad (C2)$$

$$\bar{\mathbf{\Pi}} : \frac{d\bar{\mathbf{\tau}}}{d\tau} = 2\mathbf{e}_b \cdot \bar{\mathbf{\Pi}} \cdot \mathbf{V}, \qquad (C3)$$

where we introduced the vector

$$\mathbf{V} = \hat{b}_x \frac{\partial \mathbf{u}}{\partial \hat{x}} - \frac{1}{\hat{b}} \mathbf{e}_x \times \frac{\partial \hat{\mathbf{E}}_H}{\partial \hat{x}} \qquad (C4)$$

and where $\hat{\mathbf{E}}_H$ is the electric field from the Hall and electron pressure contributions. It has been normalized with $c/(v_{0x}B_0 \sin \theta)$, and it takes the form

$$\hat{\mathbf{E}}_H = \frac{\tan \theta}{\hat{\rho}} \left( \mathbf{e}_x \times \frac{\partial \hat{\mathbf{B}}}{\partial \hat{x}} \right) \times \hat{\mathbf{B}} + \frac{M_e}{M_A \sin \theta \cos \theta} \frac{1}{\hat{\rho}} \frac{\partial \hat{\rho}}{\partial \hat{x}} \mathbf{e}_x. \qquad (C5)$$

This equation was found after assuming isotropic electron pressure and the equation of state introduced in Sec. II. However, the contribution of the last term in Eq. (C5) vanishes once it is inserted in Eq. (C4).

Finally, note that Eq. (C3) involves all the components of $\bar{\mathbf{\Pi}}$. The FLR pressure tensor is symmetric, hence it has six different components. The first row (column) is given
by Eq. (B6). The components $\Pi_{yy}$, $\Pi_{yz}$, and $\Pi_{zz}$ need to be derived from Eqs. (6)–(8). Similarly to Eq. (B6), they can be expressed as

$$ \Pi_{yy}e_x + \Pi_{yz}e_y + \Pi_{zz}e_z = \frac{1}{M_A \cos \theta} \frac{1}{\hat{b}} \frac{\partial \mathbf{u}}{\partial \hat{x}}, \quad (C6) $$

with the components of $\tilde{\mathbf{N}}$ given by

$$ N_{11} = -\hat{b}_x \hat{b}_y \hat{b}_z (P_{\perp} - 4a_p P_{\|}), \quad (C7) $$

$$ N_{12} = \frac{1}{2} \hat{b}_z P_{\perp} (1 + 3 \hat{b}_y^2), \quad (C8) $$

$$ N_{13} = \frac{1}{2} \hat{b}_y [P_{\perp} (4 + \hat{b}_x^2 - 4 \hat{b}_y^2 - \hat{b}_z^2) - 8a_p P_{\|} \hat{b}_x^2], \quad (C9) $$

$$ N_{21} = \frac{1}{4} \hat{b}_x [2P_{\perp} (\hat{b}_y^2 - \hat{b}_z^2) + 8a_p P_{\|} (\hat{b}_z^2 - \hat{b}_y^2)], \quad (C10) $$

$$ N_{22} = -\frac{1}{4} \hat{b}_y [P_{\perp} (5 + \hat{b}_x^2 - \hat{b}_y^2 - 7 \hat{b}_z^2) - 8a_p P_{\|} \hat{b}_x^2], \quad (C11) $$

$$ N_{23} = \frac{1}{4} \hat{b}_z [P_{\perp} (5 + \hat{b}_x^2 - 7 \hat{b}_y^2 - \hat{b}_z^2) - 8a_p P_{\|} \hat{b}_x^2], \quad (C12) $$

$$ N_{31} = \hat{b}_x \hat{b}_y \hat{b}_z (P_{\perp} - 4a_p P_{\|}), \quad (C13) $$

$$ N_{32} = -\frac{1}{2} \hat{b}_z [P_{\perp} (4 + \hat{b}_x^2 - \hat{b}_y^2 - 4 \hat{b}_z^2) - 8a_p P_{\|} \hat{b}_x^2], \quad (C14) $$

$$ N_{33} = -\frac{1}{2} \hat{b}_y P_{\perp} (1 + 3 \hat{b}_z^2). \quad (C15) $$

[1] K. Stasiewicz, P. K. Shukla, G. Gustafsson, S. Buchert, B. Lavraud, B. Thidé, and Z. Klos, Phys. Rev. Lett. 90, 085002 (2003).

[2] K. Stasiewicz, M. Longmore, S. Buchert, P. K. Shukla, B. Lavraud, and J. Pickett, Geophys. Res. Lett. 30, 2241 (2003).

[3] D. Perrone, O. Alexandrova, A. Mangeney, M. Maksimovic, C. Lacombe, V. Rakoto, J. C. Kasper, and D. Jovanovic, Astrophys. J. 826, 196 (2016).

[4] B. T. Tsurutani, G. S. Lakhina, A. Sen, P. Hellinger, K.-H. Glassmeier, and A. J. Mannucci, J. Geophys. Res. (Space Phys.) 123, 2458 (2018).

[5] V. Génot, E. Budnik, P. Hellinger, T. Passot, G. Belmont, P. M. Trávníček, P. L. Sulem, E. Lucek, and I. Dandouras, Ann. Geophys. 27, 601 (2009).

[6] E. A. Kuznetsov, T. Passot, and P. L. Sulem, Phys. Rev. Lett. 98, 235003 (2007).

[7] H. Kever and G. K. Morikawa, Phys. Fluids 12, 2090 (1969).

[8] T. Kakutani and H. Ono, J. Phys. Soc. Jpn. 26, 1305 (1969).

[9] A. Rogister, Phys. Fluids 14, 2733 (1971).

[10] E. Mjølhus, J. Plasma Phys. 16, 321 (1976).

[11] E. Mjølhus and J. Wyller, J. Plasma Phys. 40, 299 (1988).

[12] T. Hada, Geophys. Res. Lett. 20, 2415 (1993).

[13] E. Mjølhus, Phys. Scr. T122, 135 (2006).

[14] G. Sánchez-Arriaga, D. Laveder, T. Passot, and P. L. Sulem, Phys. Rev. E 82, 016406 (2010).

[15] A. A. Mamun, Phys. Rev. E 55, 1852 (1997).

[16] M. Tribeche, R. Amour, and P. K. Shukla, Phys. Rev. E 85, 037401 (2012).

[17] G. Williams, I. Kourakis, F. Verheest, and M. A. Hellberg, Phys. Rev. E 88, 023103 (2013).

[18] L.-N. Hau and B. U. O. Sonnerup, J. Geophys. Res. 94, 6539 (1989).

[19] L.-N. Hau and B. U. O. Sonnerup, J. Geophys. Res. 95, 18791 (1990).

[20] G. Sánchez-Arriaga, Phys. Plasmas 20, 102102 (2013).

[21] L.-N. Hau and B. U. O. Sonnerup, J. Geophys. Res. 96, 15767 (1991).

[22] J. F. McKenzie and T. B. Doyle, Phys. Plasmas 8, 4367 (2001).

[23] J. F. McKenzie and T. B. Doyle, Phys. Plasmas 9, 55 (2002).

[24] E. Mjølhus, Nonlin. Proc. Geophys. 16, 251 (2009).

[25] A. R. Champneys, Physica D 112, 158 (1998).

[26] P. L. Sulem and T. Passot, J. Plasma Phys. 81, 325810103 (2015).

[27] A. Macmahon, Phys. Fluids 8, 1840 (1965).

[28] N. Yajima, Prog. Theor. Phys. 36, 1 (1966).

[29] G. Sánchez-Arriaga and E. Siminos, J. Phys. A 50, 185501 (2017).

[30] J. Harterich, Physica D 112, 187 (1998).

[31] R. Devaney, J. Diff. Eq. 21, 431 (1976).

[32] D. Laveder, L. Marradi, T. Passot, and P. Sulem, Planet. Space Sci. 59, 556 (2011).