Finite Larmor radius influence on MHD solitary waves

MHD solitons are studied in a model where the usual Hall-MHD model is extended to include the finite Larmor radius (FLR) corrections to the pressure tensor. The resulting 4-dimensional set of differential equations is treated numerically. In this extended model, the point at infinity can be of several types. Necessary for the existence of localized solutions is that it is either a saddle-saddle, a saddle-center, or, possibly, a focus-focus. In cases of saddle-center, numerical solutions for localized travelling structures have been obtained, and compared with corresponding results from the Hall-MHD model.

Correspondence to: E. Mjølhus (einar.mjolhus@uit.no)All of the theoretical results referred to above, were obtained in models where MHD equations were extended with Hall dispersion.This kind of models will henceforth be referred to as Hall-MHD models.However, it has recently been emphasized that the Hall-MHD models in many cases lack consistency (Pokhotelov et al., 2005), and that dispersive contributions of thermal origin must be included for consistency.
In the present paper, an effort to meet this requirement is presented.To be more specific: the Hall-MHD model is extended to include the Finite Larmor Radius (FLR) correction to the pressure tensor, as derived by MacMahon (1965) and also presented in Yajima (1966).
This extension is rather challenging.While the assumption of stationary travelling wave solution in Hall-MHD models leads to a system of two coupled first order differential equations plus an algebraic equation, for which a saddle connection is to be computed, the FLR extension leads to a singular system of 5 coupled differential equations.The singularity of this system is such as to produce an algebraic equation for consistency, so that the effective dimensionality of the system is four (Hau and Sonnerup, 1991).In the Hall-MHD case, the upstream state may generically be either a center or a saddle point.A necessary condition for localized solution is that it is a saddle point.In contrast, because of the increased dimensionality, in the FLR case a larger set of possibilities occurs.The cases that allow localized solutions, include saddle-center, saddle-saddle, and, possibly focus-focus.
Due to lack of adequate numerical method, it has not been possible so far to obtain localized solutions in the two latter types of situation.However, in cases of saddle-center, this has been achieved.
The plan of paper is as follows: in Sect.2, the basic model is formulated.In Sect. 3 some of its mathematical properties are described, including in particular the treatment of the singular FLR matrix.In Sect.4, the analysis of the upstream Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.state is described.In Sect.5, some aspects of the numerical procedure are discussed.In Sect.6, some numerical results are presented.In Sect.7, a perturbation method is introduced, which can also be used in saddle-saddle and focusfocus cases.Section 8 contains a concluding discussion.

FLR-Hall-MHD model
We start with the set of equations Equation ( 1) is the equation of continuity, where ρ is the mass density and v is the ion velocity.Equation ( 2) is the equation of motion, where are the ion and electron pressure tensors, and B is the magnetic field.I ∼ is the unit (identity) tensor.Equation ( 3) is the magnetic induction equation, where the last term on the righthand side is the Hall term.The constants m i , e, c are the ion mass, the electron charge, and the speed of light respectively.
For the electron pressure, we shall throughout use an assumption of scalar pressure and isothermal equation of state Here, T e is the electron temperature, which will be assumed constant, since the electron temperature is assumed constant along magnetic field lines, and κ is the Boltzmann constant.
For the ion pressure tensor, we set where the gyrotropic part is thus, p ,⊥ are the parallel and perpendicular pressures, and e b =B/B is the unit vector along the magnetic field.The tensor P

∼
(1) i then represents the FLR correction according to MacMahon (1965) and Yajima (1966).It can be expressed as follows: i,1 i,2 where In Eq. ( 8) i =eB/m i c is the ion gyro frequency, where only one kind of ions is assumed.The notation In order to complete the modelling, equations of state for p and p ⊥ must be specified.In this paper, only the doubleadiabatic model will be considered.The double-adiabatic equations of state are based upon the assumption of neglecting the parallel heat fluxes, which is also underlying the derivation of the FLR corrections (MacMahon, 1965).It has therefore been argued that it may be inconsistent to use other pressure assumptions in combination with FLR corrections (Krauss-Varban et al., 1995).Other authors have argued to use other pressure assumptions, based, for example, on observations (for example, Hau et al., 2005;Hau and Sonnerup, 1991;Sauer et al., 2007); implying, correctly, that parallel heat fluxes are often not negligible.In this paper, it has been chosen to stay with the assumption that is consistent with the derivation of the FLR corrections.This set of equations is next specialized for stationary moving localized solutions with all state variables varying in one direction only.A coordinate system (x, y, z) is defined by an orthonormal, righthanded, constant base e x , e y , e z .
Here, e x is the wave normal direction, i.e., Moreover, the y-and z-directions are chosen such that the upstream magnetic field is Without loss of generality, the upstream velocity can be chosen as Then, v 0x defines the velocity of propagation.When v 0x is negative, the structure is moving in the positive x direction, and upstream is at x→+∞.This will be assumed, and we shall put C=−v 0x for the soliton velocity.With these assumptions and definitions, our equations can be brought to the following form: Here, F=F x e x +F y e y +F z e z has components and G=G y e y +G z e z has the components The state variables u=ue x +u y e y +u z e z , and b=b y e y +b z e z are defined by the normalizations following McKenzie and Doyle (2002).By this choice, the upstream state is Moreover, by the equation of continuity, u is related to the density by and, defining P ⊥ =p ⊥ /p 0⊥ , P =p /p 0 where p 0⊥ , p 0 are the upstream values of the pressure components, the equations of state (9) can be expressed as where In Eq. ( 16), the functions P and χ are defined as and χ 0 =χ (1, 1) is the upstream value of χ.The parameters are defined as where Finally, the tensor A ∼ can be expressed on dyadic form as and In the sequel, the tensor A ∼ will be referred to as the 1-D FLR tensor.
FLR-Hall-MHD equations on the form Eqs. ( 14), (15) for 1-D stationary moving structures have earlier been developed and studied by Hau and Sonnerup (1991) for the purpose of providing a smooth model of rotational discontinuities.In that case, however, an assumption of scalar pressure was made.The 1-D FLR tensor used by Hau and Sonnerup (1991) agrees with the present result Eqs. ( 28 14), (15) looks like a 5-dimensional system of first order differential equations.However, matters are made more complicated by the fact that the tensor A ∼ is singular, as also noted by Hau and Sonnerup (1991): where the left and right null vectors are The normalizations of L and R were chosen such that L•R=1 (note that r +r ⊥ =1).
A consequence of Eq. ( 31) is the constraint where the function H (u, b) is defined by This follows readily by dotting Eq. ( 14) from the left with L. Equation ( 33) defines a 4-dimensional surface in the 5dimensonal u, b space.This surface is the state space of Eqs. ( 14), (15).We proceed to describe how to make Eq.( 14) explicit.To this end, the eigenvectors of A ∼ are utilized.In addition to the eigenvalue 0 with left and right eigenvectors L and R, A ∼ has eigenvalues ±iµ, where µ is defined by Eq. ( 32).The vectors and these vectors satisfy the orthogonality and normalization conditions Then, Eq. ( 14) implies where has to satisfy We decompose F in the base R, S, T. However, the component of F along R must vanish because of Eqs. ( 31) and ( 36).Therefore where, by Eqs. ( 39), Inserting Eqs. ( 41) and ( 43) into Eq.( 42) and using Eqs.( 36), ( 38), (39), gives while w R must be determined from the condition of maintaining Eq. ( 33).This condition can be derived as follows: since H remains zero, one has dH /dx=0 implying which, by Eq. ( 41) gives where and noting that (∂H /∂u)•T=0.It follows from Eq. ( 46) that solutions cannot cross the set on the surface H (U )=0 in U space defined by This condition plays a role similar to the sonic circle occurring in the Hall-MHD theory (see, e.g., Mjølhus, 2006), although, the location of this set is not necessarily approximating the sonic circle.

Necessary condition at the upstream state
The upstream state u 0 =(1, 0, 0), b 0 =(0, 1) is an equilibrium point for Eq. ( 14), (15), i.e., F(u 0 , b 0 )=0, G(u 0 , b 0 )=0.A localized solution to Eqs. ( 14), ( 15) must start and end in this state as x→±∞.Therefore, a necessary condition for such solutions to exist, is that the linearization of Eqs. ( 14), ( 15) around the upstream state permits both growing and decaying solutions.In the Hall-MHD case, being two-dimensional, this amounts to a condition that the upstream state be a saddle point.The present case, being essentially 4-dimensional, allows more possibilities.Denoting linearizing around the upstream state, and assuming with where the 5×5 matrices M and N are built up as follows: where M 11 =(∂F/∂u) 0 (3×3), M 12 =(∂G/∂u) 0 (2×3), M 21 =(∂F/∂b) 0 (3×2), and M 22 =(∂G/∂b) 0 (2×2), and where A 0 is the matrix of the tensor A ∼ (referred to the base e x , e y , e z ), and I 2 is the 2×2 unit matrix.The condition for nontrivial V is det K=0, which gives the biquadratic with κ=k .The coefficients are The coefficient p 0 vanishes when the velocity C=|v 0x | is equal to one of the MHD velocities of propagation (intermediate velocity, d i =0, or slow/fast magnetosonic velocity, d ms =0).The coefficient p 2 is given by where µ 0 and R0 are the values of µ and R at the upstream state.For p 1 we have where thus, C s is the acoustic velocity of double adiabatic MHD.The coefficients q 1 and q 2 are given by complicated expressions which are given in the appendix.It can be seen that δ 0 =0 reproduces the Hall-MHD result (e.g., Mjølhus, 2006).
The classification of the upstream state is thus determined by the biquadratic P (κ) in Eq. ( 52).The biquadratic structure, which reflects the reversibility described in the next section (Symmetry 1), limits the possibilities.The generic cases are: (i).Saddle-center: This is when the two roots for κ 2 have opposite signs; then the four roots for κ are one positive, one negative and two purely imaginary complex conjugates.Consequently, the tangent space of the upstream state, represented by the linearization above, has one expanding (κ>0) and one contracting (κ<0) eigenvector direction Eq. ( 48), and a two-dimensional subspace of oscillations.
(ii).Saddle-saddle: This is when the two roots for κ 2 are both positive; then the four roots for κ are two positive and two negative.In the tangent space of the upstream state, there is a two-dimensional expanding subspace and a two-dimensional contracting subspace.
(iii).Center-center: This is when the two roots for κ 2 are both negative.Then, the four roots for κ are two pairs of complex conjugate purely imaginary roots.Then the whole tangent space of the upstream state consists of oscillations.
Referring now to Eq. ( 52), when p 0 and p 2 have opposite signs, the upstream state is a saddle-center.When they have the same sign and still p 2 1 −4p 0 p 2 >0, it is saddle-saddle when p 1 has the sign opposite of the sign of p 0 , p 2 , and center-center when p 1 has the same sign as p 0 , p 2 .Finally, when p 2 1 −4p 0 p 2 <0 it is focus-focus.The necessary condition is satisfied when we have saddlesaddle or saddle-center, while it is not satisfied when we have center-center.Cases of focus-focus have not been investigated within this work, although there are cases of this in some parameter ranges where the Hall-MHD model yields a solution.
In the Hall-MHD case, the necessary condition was determined by the velocity of propagation relative to the MHD velocities and the sonic velocity C S .The first of these depend on p 0 , which thus changes sign at the same thresholds as before.The coefficient p 1 changes sign approximately where q 0 =0, i.e., C 2 s =C 2 , when |δ 0 | is small.An example of the classification is shown in Fig. 1.In Fig. 2, the various MHD velocities have been plotted, and comparing the two figures, the significance of the various regions in Fig. 1 can be understood.In addition to the fast and slow MHD velocities (blue), the intermediate velocity (red), and the sonic velocity C S given by Eq. ( 56) (green), also the velocity C FLR at which 0R =0, has been plotted (magenta).The latter defines the locus at which the coefficient p 2 of Eq. ( 52) vanishes.
In Fig. 2, four examples of maximal existence intervals for soliton families (Mjølhus, 2006) are indicated by black color.We shall return to their interpretation in Sect.6.Here, we shall be content to note that these intervals also satisfy the necessary condition for existence of localized stationary solution in our FLR-Hall-MHD model, either exact (1,3,4) or approximately (2).

Numerical procedure
No proof has been obtained that exact solutions in the form of solitary waves exist for the 4-dimensional system Eqs.( 14 computed solutions will be presented.In the present section, some of the numerical challenges will be discussed. The integration of Eq. ( 14), ( 15) proceeded by a second order Adams-Bashforth method, starting up with one step Euler's method.The first challenge to meet, was to treat the singularity of the tensor A ∼ entering Eq. ( 14), which implies that Eq. ( 14) is not on explicit form.In the numerical procedure, Eq. ( 14) is replaced by Eq. ( 40), where w is calculated as described in Sect.3. The following shortcut was made: b entering Eq. ( 46) was calculated as b =δH /δx, where δH is the residual value of H after the previous step, and δx is the step length.
In order to numerically construct a valid localized solution, it is made use of the following reversibility property of Eq. ( 14), (15): Here, the operator J 1 is defined by A consequence of this is that if there exists a solution U 1 (x) for x<0 with u y (0)=0, b y (0)=0, which approaches the upstream state as x→−∞, this solution can be joined with U 2 (x)=(J 1 U 1 )(x) for x>0 to produce a solution which approaches the upstream state as x→+∞.Then is a localized solution for −∞<x<∞.(Here, we exploit that a solution U (x) is uniquely determined by its value at one point, U (x 0 ); in addition, this initial point must be on the constraint manifold Eq. 33.) Numerically, this can be done as follows: (i) Start a numerical solution at a point U init near the upstream state U 0 , on (or near) the unstable manifold; i.e., where τ is small and V satisfies Eq. ( 48) with a positive value of k satisfying Eq. ( 52).In practice, τ =10 −3 has been mostly used with V having unit Euclidean norm, and the O(τ 2 ) contributions to the exact unstable manifold were neglected.
(ii) Continue the computation until b y changes sign.
(iii) Denote the numerical solution generated this way Û (x), and assume that b y changes sign at x=x c .Now, if also u y changes sign at x=x c (or, numerically, very nearby), then one can use the translational invariance to define U 1 (x)= Û (x+x c ), where U 1 (x) is now defined numerically on a grid for −x c <x<0.Then U 2 (x)=(J 1 U 1 )(x) will be defined on a numerical grid for 0<x<x c , and will be a numerical approximation defined on a grid from −x c to x c , ending near U 0 at ±∞.
In the saddle-saddle cases, this procedure could be thought of as the basis for a shooting method.Namely, in those cases, there are two independent tangent vectors, V1 and V2 , spanning tangent space of the unstable manifold.By starting from a sequence of linear combinations cos ϕ V1 + sin ϕ V2 , one could think of numerically determining a value ϕ c for which u y (x c ) 0.
This procedure has, however, not been successful.In most saddle-saddle cases tried numerically, the numerical solution went to overflow.In some cases, R =0 was encountered.At low beta, this is not unreasonable: at low beta, one of the two spanning tangent directions for expanding solutions, V1 , approximates the one valid for the Hall-MHD model.The other one, V2 , must necessarily have a large value of the expansion rate, k 2 , tending to infinity as beta goes to zero (δ 0 →0 in Eq. 52).Thus, starting on the V1 direction in a low beta case, small numerical errors will make the far more rapid expansion along V2 take over.In the cases encountered, this growth did not stop, and thus went to overflow.
In the saddle-center cases, the situation is different: only two numerical solutions generated the way described above, are possible.The two correspond to two possible orientations of V , or equivalently, two possible signs of τ , referring to Eq. ( 57).The two possibilities correspond to the possibilities of generating numerically a bright vs. a dark soliton, where bright vs. dark refer to whether b 2 >1 or b 2 <1 throughout the soliton (Mjølhus, 2006).They are identified by the sign of the perturbation of b z ; if the latter is positive (negative), a bright (dark) soliton will result, provided a solution exists.Now, computing a numerical solution up to the point x c as described above in a saddle-center case, since there is no parameter to vary, it is not obvious that the condition u y (x c )=0 is hit.However, in many cases, very small values of u y (x c ) were obtained, indicating numerically that a soliton solution exists.Examples will be shown in the next section.Briefly, this happened for parameter values at which a solution exists in the Hall-MHD case, although exceptions may occur.In the next section, this will be further discussed through examples.
In some cases, we have gone a step further, in order to check this criterion of existence: if the exact value of u y (x c ) is zero, then a nonzero numerical value is due to numerical errors.These have two sources: (i) the use of the tangent to the unstable manifold instead of the exact unstable manifold for the initial condition; this error depends on the numerical parameter τ as O(τ ) 2 (cf.Eq. 59).(ii) The other source is the integration error.Since our integration method is of second order, the global integration error is expected to depend on the step length δx as (δx) 2 .
Assume that the exact value of u y (x c ) vanishes.Then, if the computation is repeated for decreasing values of δx and τ , the numerical value of |u y (x c )| is expected to decrease as the square of these parameters.On the other hand, if the exact value |u y (x c )| =0 and large enough to be numerically detectable, then the values of |u y (x c )| obtained by the above procedure will settle on a constant value even when the numerical parameters are further decreased.
The value of u y (x c ) used in this procedure is obtained by linear interpolation, using the two last values of b y , and the corresponding two last values of u y .
In practice, a decay as (δx) 2 was not achieved; instead, |u y (x c )| decayed as |δx|.A clear explanation of this has not been obtained; the most likely candidate is the shortcut by which the quantity b was calculated, as described in the beginning of this section.However, we have not succeeded in replacing this with anything better.

Numerical examples
Let us recall from Mjølhus (2006) that the soliton solutions of the Hall-MHD model come in families, characterized by velocity-amplitude relationships.For fixed wave normal direction θ, these occur within propagation velocity intervals for which the necessary condition for existence of stationary localized solutions are satisfied.Four such intervals, marked 1-4, are indicated by black vertical broken lines in Fig. 2. All of these sets also satisfy the necessary condition in the FLR-Hall-MHD model (though, interval 2 only approximately), as is easily seen by comparing Figs. 2 and 1.Within each of these maximal existence intervals, there is potentially two families of solitons, namely a dark, carrying a dip in b 2 , and a bright, carrying an increase in b 2 .In addition to the necessary condition for decaying solutions at the upstream state, a condition that the solution curve must not  Fig. 3.An example of a solution obtained by the method described in Sect. 5.The step length of integration δx=0.001, and τ =0.001 (entering Eq. 59).The residual value of u y (x c ) −4.6×10 −7 .The parameters were as indicated by the diamond in Fig. 2, with soliton velocity C=0.9 v.The length unit in this figure as well as Figs. 5 and 8, is the dispersive length .
encounter the sonic circle (Mjølhus, 2006) also had to be satisfied.Then, only a subset of these potential families actually exist as Hall-MHD soliton solutions.The maximal existence interval marked 1 in Fig. 2, contains a full Hall-MHD dark family, with a small-amplitude limit when C→C fast .It is an example of a fast magnetosonic soliton family, cold case.From Fig. 1, it can be seen that the upstream state is of saddle-center type in this case, and so we have a case in which the numerical method may work.In this case, the FLR-Hall-MHD model also produces the full fast magnetosonic (dark) family.A typical case of a solution in this family is shown in Fig. 3.The green curve was obtained by a perturbation method to be described in Sect.7.
In Fig. 4, the velocity-amplitude relationship (C versus min(b 2 )), as computed from the two models, is plotted.They are observed to coincide.This does not appear to be an exact relationship between the two models; at higher beta (v ⊥ 1), the numerical results showed discrepancy between the two.
The existence interval of Hall-MHD bright solitons in the maximal existence interval 1 of Fig. 2, is very tiny, though nonvanishing.In Mjølhus (2006), it was demonstrated that a double family of dark and bright nearly circularly polarized solitons exists sufficiently near the intermediate velocity.This double family was termed the Alfvenic soliton family.It is implied that the dark branch of this double family is the same as the magnetosonic family, as argued in Mjølhus (2006).In the Hall-MHD model, the value of C had to be extremely close to the intermediate velocity.With C int =v A cos θ=0.8660 (putting v A =1), a Hall-MHD bright soliton was numerically obtained for C=0.86625, while for C=0.8663 the solution curve encountered the sonic circle.In the present FLR-Hall-MHD model, a solution was, in fact, numerically obtained in a larger velocity interval above C int , although still tiny.Figure 5 shows the case with C=0.867.For C=0.868, the locus R =0 was encountered (cf.Sect.3). Figure 6 shows the computed velocity-amplitude relationship for this double Alfvenic family.The stars represent computations in the Hall-MHD model.Again, it is observed that the two models yield practically the same velocity-amplitude curve in terms of b 2 .In Fig. 7, the hodograms for the two solutions indicated by circles in Fig. 6, are exhibited, in order to demonstrate the circular polarization of solitons in this Alfvenic family.
The existence intervals 2 and 4 are partially of type saddlesaddle and partially of type focus-focus, and so, FLR solutions have not been obtained in these ranges.In particular, interval 4 is in a range which is potentially of physical interest, since it contains a slow magnetosonic family (warm case) (Mjølhus, 2006), which is the kind of solitons invoked in the context of magnetic holes (Baumgärtel, 1999;Stasiewicz et al., 2003;Stasiewicz, 2004).
Existence interval 3 is an existence interval bounded by C fast and the sonic velocity C s .In Hall-MHD, only a dark family can potentially exist in this case (Mjølhus, 2006), termed the fast magnetosonic family, warm case.In this case, our conclusion is that this family does not exist in the FLR-Hall-MHD model, for most of the maximal existence interval, though, the numerics indicate that it does exist for C sufficiently close to C fast .In part of the interval, the value of u y (x c ) was small, but applying the procedure described at the end of Sect.5, this value did not decrease when the numerical parameters δx, τ were decreased.However, at C=0.9378 C fast (=1.1 v A ), a residual value of u y (x c ) could not be detected when δx was halved 16 times, starting with δx=0.01.Similar behaviour has been found at other values of v ⊥ and θ.This indicates that the FLR-Hall-MHD version of the fast magnetosonic soliton (warm case) at most exists for C sufficiently close to the small-amplitude end C fast .The example C=1.1 v A is shown in Fig. 8, while a velocityamplitude relationship obtained in the range in which FLR solution appears to exist, is shown in Fig. 9.It remains to discuss cases exemplified by region 5 of Fig. 2: in this case, the upstream state is of type center-saddle according to the FLR-Hall-MHD model.However, in this case, the upstream state is of type center in the Hall-MHD model, so no solitons exist in Hall-MHD in that range.In the FLR model, computations either hit the set R =0, or a nonzero value of u y (x c ) was obtained.It is concluded that no new type of MHD soliton was found in this range in the FLR-Hall-MHD model.

Perturbation method
At low beta, it is possible to treat δ 0 , i.e., the upstream value of δ entering Eq. ( 14) and defined in Eq. ( 27), as a small parameter.Then, one can treat the FLR effects as a perturbation, starting with the Hall-MHD model as zeroth order.It was found that, in order to do so, the method of strained coordinates (Nayfeh, 1973) had to be included, because the perturbation will generally influence the width of the soliton solutions.Thus, the following expansions are assumed: and The latter describes a near-identity coordinate transformation from old coordinate x to new coordinate x 0 , and where the coefficients s 1 , s 2 , . . ., are at our disposal.This transformation implies the differentiation rule Introducing these expansions, one first gets from Eq. ( 14): In the Hall-MHD description, Eq. ( 64) is solved for u in terms of b: Inserting this into Eq.( 15) gives Equation ( 66) is the Hall-MHD model.It is a set of equations for the two transverse components of the magnetic field.Next, Eqs. ( 14) and ( 15) to order δ 0 read Here, D u F (0) , D b F (0) , D u G (0) , D b G (0) are the respective Jacobian matrices taken at the unperturbed solution; δ (0) is δ as given by Eq. ( 27), taken at the unperturbed solution, while A ∼ (0) is the tensor A ∼ given by Eqs. ( 28)-( 30) taken at the unperturbed solution.
Provided D u F (0) is invertible, Eq. ( 67) can be solved for u (1) .This can be inserted into Eq.( 68) to yield the linear, inhomogeneous differential equation for the perturbation of the magnetic field b (1) .Here, the 2×2 matrix P ∼ and the 2-vector Q are Next, the invertibility of D u F (0) is discussed.Invertibility of this matrix is lost precisely when the sonic circle [4] is encountered.Thus, when a Hall-MHD solution exists within a maximal existence interval, D u F (0) is invertible along the unperturbed solution.
The parameter s 1 determining the first order coordinate straining, is determined by the condition b (1)  y (x 0c ) = 0 (72) Here, x 0c is the value of x 0 at which component b (0) y of the Hall-MHD solution changes sign, i.e., the "half-way" point, similar to what was described in Sect. 5 for the FLR-Hall-MHD model.
In practice, solutions to the Hall-MHD model are determined numerically.The value of s 1 was determined by first solving Eq. ( 66) to determine the Hall-MHD solution b (0) (x 0 ) and the "half-way" point x 0c .Then Eq. ( 69) was solved numerically for two different values of s 1 , inserting the numerically obtained b (0) (x 0 ).This produces two values of b (1) y (x 0c ), generally violating Eq. ( 72).A secant method calculation using these two values of b (1) y (x 0c ) then produces a new value of s 1 .It can readily be demonstrated that b (1) (x 0 ) depends linearly on s 1 .Therefore, the secant calculation gives the correct value of s 1 up to numerical errors.Then, the correct b (1) (x 0 ) was determined by solving Eq. ( 69) with this value of s 1 .
In Figs. 3 and 8, the green curve shows the solution computed this way.One can observe that for the parameter sets of these cases, the solution calculated by the perturbation method practically overlaps the solution calculated from the full FLR-Hall-MHD model.In the case exhibited in Fig. 5, the solution calculated by the perturbation method deviated substantially from the FLR solution, and it is not shown here.

Discussion
In this paper, stationary localized travelling wave solutions (solitary waves, or "solitons") for models representing dispersive extensions of MHD models, have been discussed.Focus has been on extending the Hall-MHD model (McKenzie and Doyle, 2002) with dispersive effects of thermal origin.Such effects are embodied in the finite Larmor radius extension to the pressure tensor (MacMahon, 1965;Yajima, 1966).The relevant set of ordinary differential equations including this has been developed.The upstream state is an equilibrium point for this set of equations.Localized solutions are saddle connections for the upstream state.This implies that the upstream state must be a saddle-saddle, focusfocus, or a saddle-center.Numerical solutions have been obtained in saddle-center cases for which a saddle connection exists also in the Hall-MHD case (or for parameter sets nearby, cf. the case of Fig. 5), while in saddle-center cases where a saddle connection does not exist in the Hall-MHD model, the numerics showed that it does not exist in the FLR-Hall-MHD model either.In saddle-saddle cases, a numerically constructed saddle connection was not achieved.There is presently no reason to interpret this as nonexistence of solution; it may equally well be because the numerical approach so far attempted, is inadequate.
Neither should the numerical solutions that could be obtained in the saddle-center cases be taken as a proof that a solution exists.Only a mathematical proof can decide that.Unfortunately, this has not been within reach so far.
Since the set of differential Eqs. ( 14), (15) has been obtained from a model containing no irreversible effects, one might expect it to have some Hamiltonian structure.Such a structure could have been useful in the search for a precise mathematical proof of existence of saddle connections.However, no such structure has so far been found.[A referee has made the following remark: even an energy conservation theorem is not known for this model.Such a theorem can be obtained (i) for the Hall-MHD model with isotropic pressure, (ii) MHD with double-adiabatic pressure tensor, but not for Hall-MHD with double-adiabatic pressure, so it becomes problematic even before the inclusion of FLR.This may be connected with the fact that the double-adiabatic closure builds on the frozen-in condition, which is broken when the Hall term is included.Obviously, the closure of highermoment models still deserves some attention.] One of the motivations of the present work was to throw light on the proposed theoretical explanation of magnetic holes in terms of propagating slow magnetosonic dark solitons (Baumgärtel, 1999;Stasiewicz et al., 2003).In particular, there was a need to meet the challenge of Pokhotelov et al. (2005) to study such solutions including dispersive terms due to finite ion temperature.Although this challenge has been partially met by the present contribution, there is no reason to claim that it throws any new light on the explanation of magnetic holes.Even so, some remarks to this topic are next given: it should first be mentioned that in much of the current discussion, the observed magnetic holes are associated with the mirror instability (Chandrasekhar et al., 1958;Pokhotelov et al., 2004), because they are predominantly observed in plasma states near the threshold for this instability (e.g., Sperveslage et al., 2000).Therefore, they are also sometimes referred to as mirror structures.
Since dark solitons of the one-parameter type (Mjølhus, 2006) discussed in this work have not usually been associated with instabilities, their role in this context should at most be understood as aftermaths after the instability has played itself out, and the plasma containing the mirror structures has reached a stable state.This has been the view of the present author.The question then becomes whether theory can provide stationary structures as outcome of the mirror instability, and which are clearly distinct from the slow magnetosonic solitons that have been proposed (Baumgärtel, 1999;Stasiewicz et al., 2003).Both in the recent discussion (Stasiewicz, 2005) as well as elsewhere (Treumann et al., 2004), it was claimed that such an alternative does not exist.
In the present author's view, this situation has now changed.Recently, localized static (i.e., nonpropagating) structures were theoretically obtained in fully time dependent numerical solutions to a dispersive-MHD model including both finite Larmor radius dispersion and Landau damping (Borgogno et al., 2007).The numerical run was for a set of parameter values where the pressure anisotropy p ⊥0 /p 0 was slightly above threshold for the mirror instability.In a particular case where the boundary conditions were chosen so as to maintain the mirror instability, a nonchanging and nonpropagating structure developed from initial noise (e.g., Fig. 4 of Borgogno et al., 2007) having the form of magnetic holes, anticorrelated with the density perturbations.
where b = bx e x b⊥ = by e y + bz e z b = b + b⊥ by,z = b y,z sin θ/ b = B y,z /B bx = cos θ/ b = B x /B (29) Fig. 1.Example of the classification of the upstream state.The electron temperature satisfies v te =v ⊥ .

Fig. 2 .
Fig. 2. Example of the various critical velocities entering the classification of the upstream state.The parameters are the same as for Fig. 1.The values of the velocities are in units of v A .

Fig. 4 .
Fig. 4. Velocity-amplitude relationship for the fast magnetosonic family in the case represented by 1 in Fig. 2. The stars represent results of the Hall-MHD model.

Fig. 5 .Fig. 6 .
Fig.5.An example of a bright soliton, obtained close to C int , the lower end of the existence interval 1 of Fig.2.The numerical parameters were as in Fig.3; the length unit likewise.The residual value was u y (x c ) −2.8×10 −8 .

Fig. 8 .
Fig.8.Solutions obtained by FLR-Hall-MHD as well as Hall-MHD and the perturbation method of Sect.7, for a case in the fast magnetosonic family, warm case.

Fig. 9 .
Fig. 9. Velocity-amplitude relationship for the part of the range 3 (cf.Fig. 2) in which the numerics indicate that a solution exists in FLR-Hall-MHD.