Contents

• Hemispherical transport equation: modeling of quasiparallel collisionless shocks
• 1 Introduction
• 2 Basic equations
• 3 Self-consistent closure
• 4 Application for acceleration at the plane shock
• 5 Numeric modeling of collisionless shocks
• 6 Conclusion
• References
• 7 Online Material
• Appendix A: Formal inclusion of Alfvén heating
• Appendix B: Numeric method

A&A 466, 1-9 (2007)
DOI: 10.1051/0004-6361:20065714

Hemispherical transport equation: modeling of quasiparallel collisionless shocks

V. N. Zirakashvili^1,2

1 - Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radiowave Propagation, 142190 Troitsk, Moscow Region, Russia
2 - Max-Planck-Institut für Kernphysik, 69029 Heidelberg, Postfach 103980, Germany

Received 29 May 2006 / Accepted 2 November 2006

Abstract
Using a so-called hemispherical model we derive a general transport equation for cosmic ray and thermal particles scattered in pitch angle by magnetic inhomogeneities in a moving collisionless plasma. The weak scattering through 90 degrees results in isotropic particle distributions in each hemisphere. The consideration is not limited by small anisotropies and by the condition that particle velocities are higher than characteristic flow velocity differences. For high velocities and small anisotropies the standard cosmic ray transport equation is recovered. We apply the equations derived for investigation of injection and acceleration of particles by collisionless shocks.

Key words: acceleration of particles - cosmic rays

1 Introduction

Progress in cosmic ray astrophysics has been made with the introduction of the transport equation for cosmic rays (Krymsky 1964; Parker 1965; Dolginov & Toptygin 1966; Gleeson & Axford 1969). It was successfully applied to many astrophysical problems: cosmic ray modulation in the solar wind, acceleration of particles at shock fronts, cosmic ray propagation in the Galaxy etc. This equation implies that in the strong scattering approximation the particle distribution is almost isotropic, which in particular implies that velocities of particles are higher than the characteristic flow velocity differences in the plasma. However, in several circumstances these conditions are violated and more general kinetic equations should be used which also deal with the evolution of the pitch angle distribution. Since the solution of such equations is not a simple task, we would like find a description that comes close to that implied by the standard transport equation.

Energetic particles are scattered in pitch angle relative to the mean magnetic field direction by random magnetic inhomogeneities. It is well known that in the case of particle transport along the mean magnetic field a so-called problem of scattering through the pitch-angle of 90 degrees exist. According to quasilinear theory the resonant scattering of particles is weak at small pitch-angle cosines (Jokipii 1966). The main reason for this phenomenon is the small energy density of short waves propagating along the magnetic field and scattering the particles with pitch angles close to 90 degrees. These waves may be easily damped by thermal ions. Even without such a damping, the scattering in the vicinity of 90 degrees is depressed in the case of a power-law spectrum of waves that is proportional to k^-2 or steeper. Here k is the wavenumber. This problem can be avoided if one takes into account the scattering by oblique waves or nonlinear interactions, for example magnetic mirroring (see e.g. Völk 1973, 1975; Achterberg 1981). One can expect that the scattering efficiency is small near the pitch angle of 90 degrees.

Observations of suprathermal pick-up particles in the solar wind suggest that this is the case (Fisk et al. 1997). If the scattering through 90 degrees is weak, the particle distribution is almost isotropic in each hemisphere of the velocity directions parallel or antiparallel to the mean field. So, in this case the angular dependence of particle distribution is reduced to the two number densities of forward and backward moving particles. The corresponding equations were derived by Isenberg (1997) for the case of pick-up ions in the solar wind.

We use this approach and derive general transport equations for arbitrary nonrelativistic flow velocities. We use the particle distributions in the frame moving with the medium. This permits us to consider the case of arbitrary particle velocities. Since the number density of particles in both hemispheres may substantially differ, large anisotropies can also be taken into account. The equations derived have a broad range of applicability: propagation of solar energetic particles and pick-up ions in the solar wind, injection into diffusive shock acceleration and other processes.

As an example of the application of the equations derived, we consider the problem of diffusive shock acceleration. Since our equations describe propagation of thermal as well as energetic particles in the vicinity of the collisionless shock front, they can be used for the determination of the shock velocity profile and the spectra of thermal and nonthermal particles without any assumptions regarding the injection efficiency. One only has to prescribe the scattering law. In this sense the approach is similar to the one used by Ellison & Eichler (1984). They applied a Monte-Carlo technique for the solution of the kinetic equation.

The basis of such an approach is the idea that the interaction of upstream and downstream plasmas may result in instabilities similar to firehose instability (Parker 1961; Quest 1988). The magnetic inhomogeneities produced by such instabilities may play the role of scattering centers and provide the shock dissipation and heating. These ideas are also the basis of the so-called hybrid simulations of collisionless shocks (Leroy & Winske 1983; Quest 1986). The ions move in self-consistent electromagnetic fields and the electrons are considered as a charged neutralizing fluid in these simulations.

The equations are derived in the next two sections. Application to the case of collisionless shock is considered in Sects. 4 and 5. Section 6 contains a discussion of results obtained and conclusions.

2 Basic equations

We start with the kinetic equation for the momentum distribution  $f({\vec r},{\vec p})$ of charged particles

$$% \frac {\partial f}{\partial t}+{\vec v}\nabla f+ {\vec F}_L\frac {\partial f}{\partial {\vec p}}=\hat {St} f .$$ (1)

Here ${\vec v}$ and ${\vec p}$ are the particle velocity and momentum respectively, ${\vec F}_L={q\vec E}+(q/c)[{\vec v} \times {\vec B}]$ is the Lorence force and the operator $\hat {St}$ describes the pitch angle scattering of particles by random magnetic inhomogeneities in the fluid frame. We shall further assume that the mean electric field $\vec E$ can be described as

$$% {\vec E}=E_\parallel {\vec b}-\frac 1c[{\vec u}\times {\vec B}],$$ (2)

where $\vec b$ is the unit vector along the direction of the mean magnetic field $\vec B$, $\vec u$is the mass velocity and  $E_\parallel$ is the component of the mean electric field which is parallel to the mean magnetic filed.

It is convenient to perform a change of variable ${\vec p}$ to ${\vec p}'$ that formally coincides with the nonrelativistic Lorence transform

$$% {\vec p}'={\vec p}-\frac {{\vec u}}{v}p.$$ (3)

Averaging the particle gyrophase angle and neglecting terms of the order $u/c\ll 1$we find (cf. Skilling 1971; Kulsrud 1983):

$$\begin{eqnarray*}\frac {\partial f}{\partial t}+(v'\mu '{\vec b}+{\vec u})\nabla f=\hat {St} f\\ [1mm] \end{eqnarray*}$$

$$\begin{eqnarray*}+\frac {\partial f}{\partial p'}\left[ \frac {3\mu '^2-1}{2}p' ... ...1-\mu '^2}{2}p'\nabla {\vec u}-\mu 'F_\parallel \right]\\ [2mm] \end{eqnarray*}$$

$$% +\frac {\partial f}{\partial \mu '}(1-\mu '^2)\left[ \frac ... ...rac {\mu '}{2}\nabla {\vec u}-\frac {F_\parallel }{p'} \right]$$ (4)

here $\mu '$ is the cosine of the pitch angle of the particle and  $F_\parallel$ is the sum of the parallel components of the electric and inertia force:

$$% F_\parallel =qE_\parallel - {\vec b}\frac {p'}{v'}\left( \f... ...rtial {\vec u}}{\partial t}+ ({\vec u}\nabla ){\vec u}\right).$$ (5)

The scattering operator can be written as

$$% \hat {St} f=\frac {\partial }{\partial \mu '}\frac {1-\mu '^2}{2}\nu (p',\mu ') \frac {\partial f}{\partial \mu '}$$ (6)

here $\nu$ is the scattering frequency.

Let us assume that the scattering frequency is large for pitch angle cosines $\vert\mu '\vert>\mu _0$. The momentum distribution f is then isotropic in these regions of the angle space. We introduce distributions of particles in the backward and forward hemispheres $N_\pm$: $f(p',\mu ')=N_-(p')$ for $\mu '<-\mu _0$ and $f(p',\mu ')=N_+(p')$ for $\mu '>\mu _0$. Integration of Eq. (4) from -1 to $-\mu _0$ and from $\mu _0$ to 1 gives the equations for these distributions:

$$\begin{eqnarray*}\frac {\partial N_\pm }{\partial t}+ \left( {\vec u}\pm \frac {... ...ec \nabla {\vec u}}\right) \frac {\partial N_\pm }{\partial p'}= \end{eqnarray*}$$

$$% \mp \left. \frac {\nu }{2}\frac {\partial f}{\partial \mu '}\right\vert _{\mu '=\pm \mu _0}.$$ (7)

It was assumed that $\mu _0 \ll 1$. Now we calculate the right-hand side of these equations. For this we must find the solution of Eq. (4) in the region $\vert\mu '\vert<\mu _0$. Since $\mu _0 \ll 1$, we can leave only the terms with derivatives on $\mu '$ in this equation with the result:

$$% -\left( \frac {v'}{2}\nabla {\vec b}+\frac {F_\parallel }{p... ...l \mu '}\frac {\nu}{2} \frac {\partial f}{\partial \mu '}\cdot$$ (8)

This equation should be solved with the boundary conditions $f(p',\pm \mu _0)=N_\pm (p')$. Substitution of the solution into the Eq. (7) then gives

$$\begin{eqnarray*}\frac {\partial N_\pm }{\partial t}+ \left( {\vec u}\pm \frac {... ...\bf\nabla} {\vec u}\right) \frac {\partial N_\pm }{\partial p'}= \end{eqnarray*}$$

$$% \mp \nu _\mp (N_+-N_-).$$ (9)

Here the frequencies $\nu _\pm$ describe the rate of particle exchange between forward and backward hemispheres:

$$% \nu _\pm (p')=\pm \frac {\frac {F_\parallel }{p'}+\frac {v'... ...frac {2F_\parallel }{v'p'}+\nabla {\vec b}\right) \right] -1},$$ (10)

and the mean free path $\lambda$ of the particle in the uniform medium is given by the expression:

$$% \lambda =v'\int ^{\mu _0}_{-\mu _0}\frac {{\rm d}\mu '}{\nu (p',\mu ')}\cdot$$ (11)

A very similar equation was derived by Isenberg (1997) for pick-up ions in the solar wind. We found more accurate expressions for the frequencies $\nu _\pm$and these coincide with the result of Kota (2000) obtained for the case $F_\parallel =0$.

We neglected the derivatives in time and space in Eq. (8). This assumption is valid if the characteristic time $\tau$ and the characteristic scale l of the problem considered are large enough: $\tau , \ l/v' \gg \mu _0^2/\nu \sim \mu _0\lambda /v'$. Since $\mu _0 \ll 1$, our derivation is valid even if the scale l is comparable with the mean free path $\lambda$. Equation (9) is exact in the mathematical limit $\mu _0\to 0$, $\lambda \to {\rm const}$.

Equations (9) can be rewritten in the conservative form:

$$\begin{eqnarray*}\frac {\partial N_\pm }{\partial t}+ \nabla \left( {\vec u}\pm \frac {v'}{2}{\vec b}\right) N_\pm % +\frac {1}{p'^2}\frac {\partial }{\partial p'}p'^2 \left( \... ...nabla} {\vec u}\right) N_\pm = \nu _\mp N_\mp -\nu _\pm N_\pm.$$ (12)

Equations derived have a simple physical meaning. The second term in the left hand side describes the transport of the particles by the medium moving with velocity $\vec u$ and the proper motion of particles along the magnetic field with the speed $\pm v'/2$ (to be compared with "coherent'' propagation of solar energetic particles considered by Earl 1974). The third term describes the energy losses or gains in the inhomogeneous flow and electric field. The right hand side corresponds to the exchange of particles between forward and backward hemispheres.

It is clear that the system of Eq. (12) is hyperbolic. It can be reduced to uncoupled equations for N₊ and N_-:

$$\begin{eqnarray*}\frac {\partial N_\pm }{\partial t}+{\vec u}\nabla N_\pm- \frac... ...+ \nabla \left( {\vec u}\mp \frac {v'}{2}{\vec b}\right) \right. \end{eqnarray*}$$

$$\begin{eqnarray*}\left. + \frac {1}{p'^2}\frac {\partial }{\partial p'}p'^2 \lef... ... {p'}{3}{\bf\nabla} {\vec u}\right) \right] \frac {1}{2\nu _\mp} \end{eqnarray*}$$

$$% \times\left[ \frac {\partial N_\pm }{\partial t}+ \left( {\... ...vec u}\right) \frac {\partial N_\pm }{\partial p'}\right]\cdot$$ (13)

In the case of high scattering frequencies the particle distributions in the different hemispheres are almost equal to each other, $N_+\approx N_-\approx N$. Assuming also a slow time evolution and particle velocities much larger than the medium velocity we come to the standard transport equation for cosmic rays (Krymsky 1964; Parker 1965; Dolginov & Toptygin 1966; Gleeson & Axford 1969):

$$% \frac {\partial N}{\partial t}+{\vec u}\nabla N- \frac {p'}... ...la} {\vec u}= (\nabla {\vec b})D_\parallel ({\vec b}\nabla )N.$$ (14)

Here the diffusion coefficient along the magnetic field $D_\parallel =v'\lambda /4$.

Another interesting case is when the energy changes and transport by the medium are negligible. Then Eq. (13) is reduced to the so-called telegraph equation (cf. Fisk & Axford 1969):

$$% \frac {\partial N_\pm }{\partial t}= -\left[ \frac {\partia... ...}{\partial t} \pm \frac {v'}{2}({\vec b}\nabla )N_\pm \right].$$ (15)

This equation describes propagation of cosmic ray particles with finite velocity and can be reduced to the diffusion equation only in the case of a slow time evolution.

3 Self-consistent closure

The equations derived in the previous section contain prescribed values for the flow velocities $\vec u$, the parallel force  $F_\parallel$ and the magnetic field $\vec B$. Such a description is a good approximation for the investigation of propagation of test particles. However, in many interesting cases these quantities should be determined self-consistently, involving the distributions N₊ and N_-.

We treat electrons as a massless fluid. In this case the Euler equation for electrons is reduced to an expression for the electric field:

$$% {\vec E +\delta {\vec E}}= -\frac {\nabla P_{\rm e}}{qn_{\r... ...c 1c[{\vec u}_{\rm e}\times ({\vec B} +{\bf\delta} {\vec B})],$$ (16)

where ${\vec u}_{\rm e}$ and $P_{\rm e}$ are the velocity and the pressure of the electron fluid respectively, $n_{\rm e}$ is the electron number density, ${\bf\delta} {\vec E}$ and  $\bf\delta {\vec B}$ are the random components of the electric and the magnetic fields respectively. The velocity  ${\vec u}_{\rm e}$ can be found from Maxwell's equation

$$% [{\bf\nabla} \times ({\vec B}+{\bf\delta} {\vec B})]=\frac ... ... q}{c} \int {\rm d}^3p({\vec v}-{\vec u}_{\rm e})(f+\delta f),$$ (17)

where $\delta f$ is the random component of the ion distribution, and we have neglected the displacement current by assuming a slow time evolution. It was also assumed that the plasma is pure hydrogen and quasineutral. Substitution of  ${\vec u}_{\rm e}$ from this equation into Eq. (16) and averaging gives:

$$\begin{eqnarray*}{\vec E}= -\frac 1c[{\vec u}\times {\vec B}]- \frac {1}{qn_{\rm... ...vec B}\times [\nabla \times {\vec B}]]+ \nabla P_{\rm e} \right. \end{eqnarray*}$$

$$% \left. +\int {\rm d}^3p{\vec p}\hat {St}f \right).$$ (18)

Here $n_{\rm i}$ is the mean number density of ions and we have used the formal definition of the scattering operator $\hat {St}f=-\frac qc \langle \frac {\partial }{\partial {\vec p}}[{(\vec v}-{\vec u}) \times {\delta \vec B}]\delta f \rangle$, where the angle brackets mean the average over the random fluctuations of  $\delta \vec B$ and $\delta f$ of the magnetic field and the momentum distribution. We also neglected the magnetic tension of the random component. Generally speaking, the scattering operator does not conserve momentum. Since the magnetic inhomogeneities are frozen into the electron fluid, this additional force acts on the electron fluid. As a result, an additional electric field due to the charge separation appears. It is described by the last term in the parentheses of Eq. (18). Such a field is indeed observed in hybrid simulations of collisionless shocks (Quest 1988).

We multiply Eq. (1) containing the electric field (18) by  ${\vec p}{\rm d}^3p$ and integrate over momentum space. This gives the Euler equation of motion:

$$% \rho \left( \frac {\partial {\vec u}}{\partial t}+ ({\vec u... ...\times [{\bf\nabla} \times {\vec B}]]- \nabla (P_{\rm e}+P_i).$$ (19)

Here

$$% P_i=\frac 13\int {\rm d}^3p'v'p'(N_++N_-)$$ (20)

is the pressure of ions and $\rho$ is the mean density. The parallel electric field can be found from Eq. (18).

Since the Eq. (1) with the electric field (18) conserves momentum, the same is true for Eq. (9). Let us multiply them by p', integrate over  $p'^2{\rm d}p'$ and subtract the first from the second. This gives after some algebra

$$% F_\parallel \frac {1}{n_i}\left[ \frac 23\int {\rm d}^3p'p'(\nu _++\nu _-)(N_+-N_-)+ {\vec b} \nabla P_i \right.$$ =

$$\left.+ {\vec b}\left( \frac {\partial {\vec u}}{\partial t}+ ({\... ...t {\rm d}^3p'' \left( \frac {p''}{v''}-\frac {p'}{v'}\right) (N_++N_-) \right].$$ (21)

The last term term in the right hand side of this expression can be neglected, since the force  $F_\parallel$ is essential only for nonrelativistic particles which determine the plasma density.

We shall assume an adiabatic evolution of the electron pressure

$$% \frac {\partial P_{\rm e}}{\partial t}+{\vec u}\nabla P_{\rm e}+ \frac 53P_{\rm e}\nabla {\vec u}=0,$$ (22)

and frozen-in magnetic field. The last term in Eq. (18) is essential for momentum conservation but can be neglected in Faraday's induction equation:

$$% \frac {\partial \vec B}{\partial t}= [\nabla \times [{\vec u}\times {\vec B}]].$$ (23)

We have thus obtained the closed system of Eqs. (12), (19), (22), (23) with expressions (10) for the frequencies $\nu _\pm$and Eq. (21) for the parallel force  $F_\parallel$. In the next sections we apply these equations to the combined problem of injection and nonlinear acceleration at a plane shock.

4 Application for acceleration at the plane shock

Let us consider the acceleration at the parallel one-dimensional shock. We can write the steady-state version of Eq. (9) for nonrelativistic particles in the shock frame as

$$% \left( u\pm \frac {v'}{2}\right) \frac {\partial N_\pm }{\parti... ...artial z} \right) \frac {\partial N_\pm }{\partial v'}= \mp \nu _\mp (N_+-N_-).$$ (24)

Here m is the mass of particles and the parallel force $F_\parallel =qE_\parallel -mu\partial u/\partial z$.

The characteristics of this system of equations for the case  $E_\parallel =0$are given by

$$% \left( \frac {v'}{u}\pm \frac 43+\frac {\sqrt{7}}{3}\right)... ...}}{3}\right) ^{\frac 12\pm \frac {1}{\sqrt{7}}} = {\rm const}.$$ (25)

They are shown in Fig. 1. An upstream particle begins its motion with low velocity from the upper left part of this u-v' plane and goes down the characteristics. The forward moving particles gain energy, backward moving ones lose energy. At any moment the particle may change hemisphere and will continue the motion along another set of characteristics. If the compression ratio of the shock is high enough, the forward moving particle can change hemisphere in the vicinity of the line u=v'/2 (shown by the dash-dotted line) and return upstream along the backward characteristic with energy gain. It can again change the hemisphere upstream and again move in a downstream direction etc. This may be repeated many times and the accelerated particle goes to the right beyond the part of u-v' plane shown in Fig. 1.

There is an another possibility for the initial acceleration. If the initial speed of the backward moving particle upstream is high enough, say 0.7u₁, it can reach the line u=v'/2 directly and return upstream with an energy gain. However, the initial speed is rather high. The number of such particles is small for high Mach number shocks.

Figure 1: Characteristics of the hyperbolic system of Eq. (24) for the case $E_\parallel =0$. Characteristics for forward moving particles and backward moving particles are shown by solid and dotted lines respectively.

The characteristics of Eq. (24) depend on $F_\parallel$. The characteristics for the case $F_\parallel =0$ are shown in Fig. 2. This case is close to the real situation because the two first terms in Eq. (21) almost cancel each other (see also the numeric modeling below). In other words, the electric force  $qE_\parallel$ almost compensate the inertia force $-mu\partial u/\partial z$ in the second term on the left-hand side of Eq. (24). This is not a simple coincidence. The electric force is the consequence of the momentum transfer and is approximately $-\partial P/\partial z/\rho$ that is exactly $u\partial u/\partial z$. The characteristics are determined by expressions

$$% v'^4\pm \frac 83uv'^3= {\rm const}.$$ (26)

Figure 2: Characteristics of the hyperbolic system of Eq. (24) for the case $F_\parallel =0$. Characteristics for forward moving particles and backward moving particles are shown by solid and dotted lines respectively.

We find the analytical solution of Eq. (24) for the case when the velocity profile is discontinuous: u(z)=u₁, z<0 and u(z)=u₂, z>0. In the regions of constant velocity u the system (24) can be rewritten as (to be compared with Gombosi et al. 1993)

$$u\frac {\partial N_+}{\partial z}=\frac {\partial }{\partial z} \left( {v'^2}-4u^2\right) \frac {\lambda }{4v'}\frac {\partial N_+}{\partial z}$$ (27)

$$N_-=N_++\left( 2{u}+v'\right) \frac {\lambda }{v'}\frac {\partial N_+}{\partial z}\cdot$$ (28)

It was also assumed that $E_\parallel =0$. The solution of these equations upstream (z<0) is

$$% N_\pm =N_\infty +(N_{1\pm }-N_\infty) \exp \frac {4u_1v'z/\lambda }{{v'^2}-4{u_1^2}}\cdot$$ (29)

Here $N_\infty$ is the distribution at $z=-\infty$, N₁₊ and N_1- are the distributions N₊ and N_-just upstream the shock. According to Eq. (28) they are related as

$$(v'+2u_1)(N_{1+}-N_\infty )=(v'-2u_1)(N_{1-}-N_\infty ), \ v'>2u_1$$

$$N_{1+}=N_{1-}=N_\infty ,\ v'<2u_1.$$ (30)

The solutions downstream are given by

$$% N_\pm N_{2+}\left( \frac12+\frac {v'}{4u_2}\right)+ N_{2-}\left( \frac12-\frac {v'}{4u_2}\right)$$ =

$$+\left( \frac {v'}{4u_2}\mp \frac12\right) (N_{2-}-N_{2+}) \exp \frac {4u_2v'z/\lambda }{{v'^2}-4{u_2^2}}\cdot$$ (31)

Here N₂₊ and N_2- are the distributions just downstream. They are related as

$$% N_{2+}=N_{2-}, \ v'>2u_2.$$ (32)

This means that distributions N₊ and N_- do not depend on z for v'>2u₂.

We find the relation between the upstream and downstream distributions $N_{1\pm}$ and $N_{2\pm}$. It depends on the solution of Eq. (24) in the transition region. We consider the case $F_\parallel =0$. Then the distributions upstream and downstream are related by characteristics (26):

$$% N_{1\pm }=F_\pm \left(v'^4\pm \frac 83u_1v'^3\right), \ N_{2\pm }=F_\pm \left(v'^4\pm \frac 83u_2v'^3\right).$$ (33)

Here $F_\pm$ are two functions to be determined. Thus we have the six Eqs. (30), (32), (33) for the six unknown functions $N_{1\pm}$, $N_{2\pm}$ and $F_\pm$.

Let us consider the case $N_\infty =\delta (v'-v_0)/(2v_0^2)$. Since the change of velocity in the transition region is governed by the characteristics (26), the solution is the following sum of $\delta$-functions:

$$% N_{1\pm }=A_{0}\delta (v'-v_{0})/(v_{0})^2 + \sum ^2_{j=1}\... ...nfty _{i=1}A^j_{i\pm }\delta (v'-v^j_{i})/(v^j_{i})^2,\\ [3mm]$$ (34)

$$% N_{2\pm }=\sum ^2_{j=1}\sum ^\infty _{i=0}B^j_{i\pm }\delta (v'-w^j_{i})/(w^j_{i})^2.$$ (35)

Here $A^j_{i\pm }$, $B^j_{i\pm }$, v^j_i and w^j_i are sequences to be determined. The index j may have the values 1 and 2. These values correspond to the two possibilities of injection of forward and backward moving particles. They were described after Eq. (25). The different values of the index i correspond to consequent states of the one particle moving between upstream and downstream regions of the shock.

If the upstream sequences $A^j_{i\pm }$ are known for v'>2u₁the downstream ones can be found from characteristics (26). Using the properties of $\delta$-functions we find that they are given by

$$% B^j_{i\pm }= \frac {v^j_i+2u_1}{w^j_i+2u_2}A_{i+}, \ w^j_i>2u_2$$ (36)

and the downstream velocity w^j_i is the solution of equation

$$% (v^j_i)^4+ \frac 83u_1(v^j_i)^3=(w^j_i)^4+\frac 83u_2(w^j_i)^3.$$ (37)

The next $A^j_{i+1\pm }$ are

$$% A^j_{i+1-}= \frac {w^j_i-2u_2}{v^j_{i+1}-2u_1}B^j_{i-}, \ A^j_{i+1+}=B^j_{i-}\frac {w^j_{i}-2u_2}{v^j_{i+1}+2u_1}$$ (38)

where the velocity v^j_i+1 is the solution of the equation

$$% (v^j_{i+1})^4- \frac 83u_1(v^j_{i+1})^3=(w^j_i)^4-\frac 83u_2(w^j_i)^3.$$ (39)

Equations (36)-(39) can be used recurrently to calculate the sequences  $A^j_{i\pm }$, $B^j_{i\pm }$, v^j_i and w^j_i.

Let us introduce the critical velocities $v_{c\pm}$. They are the solutions of equations:

$$% v^3_{{\rm c}\pm }(3v_{{\rm c}\pm }\pm 8u_1)=(48\pm 64)u^4_2.$$ (40)

The initial velocity v₀ and A₀=1/2 are given. There are four cases depending on the initial velocity v₀ and velocities  $v_{{\rm c}\pm}$.

1) v₀<v_c-. Forward and backward moving particles reach the downstream region with velocities smaller than 2u₂ and are further advected. The injection in acceleration does not occur. The downstream coefficients  $B^j_{0\pm}$ are given by

$$\begin{eqnarray*}B^1_{0+}=A_0\frac {v_0+2u_1}{w^1_0+2u_2}, \ B^2_{0-}=0 \end{eqnarray*}$$

$$% B^2_{0-}=A_0\frac {v_0-2u_1}{w^2_0-2u_2}, \ B^1_{0-}=0$$ (41)

where the velocities w¹₀ and w²₀ are the solutions of equations:

$$\begin{eqnarray*}(v_{0})^4+\frac 83u_1(v_{0})^3=(w^1_0)^4+\frac 83u_2(w^1_0)^3 \end{eqnarray*}$$

$$% (v_{0})^4-\frac 83u_1(v_{0})^3=(w^2_0)^4-\frac 83u_2(w^2_0)^3.$$ (42)

All other coefficients $A^j_{i\pm }=0$, $B^j_{i\pm }=0$ for i>0.

2) $v_{\rm c-}<v_0<v_{\rm c+}$. The most important case for injection. The forward moving particle reaches the downstream region with a velocity less than 2u₂ and again is not injected into acceleration. The coefficients  $B^1_{i\pm}$ and  $A^1_{i\pm}$ are the same as in the previous case.

The backward moving particle goes along characteristics and returns upstream (see Fig. 2). Using properties of $\delta$-functions in Eq. (33) and conditions (30) we found

$$% A^2_{1-}=A_{0}\left\vert \frac {v_0-2u_1}{v^2_1-2u_1}\right\vert, \ A^2_{1+}=A^2_{1-}\frac {v^2_1-2u_1}{v^2_1+2u_1}\cdot$$ (43)

Here the velocity v²₁ is the solution of the equation

$$% (v_0)^4- \frac 83u_1(v_0)^3=(v^2_1)^4-\frac 83u_1(v^2_1)^3.$$ (44)

It is clear that v²₁ is close to 8u₁/3 (see Fig. 2). The other coefficients  $A^2_{i\pm}$ and  $B^2_{i\pm}$ can be found consecutively using Eqs. (36)-(39).

3) $v_{\rm c+}<v_0<2u_1$. Forward and backward moving particles reach the downstream region with velocities greater than 2u₂ and are injected into acceleration. The coefficients  $A^2_{i\pm}$ and  $B^2_{i\pm}$ are the same as in the previous case. The coefficients  $B^1_{0\pm}$ are given by

$$% B^1_{0\pm }= \frac {v_0+2u_1}{w^1_0+2u_2}A_{0},$$ (45)

and the downstream velocity w¹₀ is the solution of the equation

$$% (v_0)^4+ \frac 83u_1(v_0)^3=(w^1_0)^4+\frac 83u_2(w^1_0)^3.$$ (46)

The coefficients $A^1_{1\pm}$ can be found from Eqs. (38) and (39). Now the coefficients  $B^1_{i\pm}$ at i>0 and  $A^1_{i\pm}$ at i>1 can be found consecutively using Eqs. (36)-(39).

4) v₀>2u₁. The coefficients $A^j_{i\pm }$ and  $B^j_{i\pm }$ at i>0can be found consecutively using Eqs. (36)-(39).

Using Eqs. (36) and (38) one can relate A^j_i+1+ and A^j_i+:

$$% A^j_{i+1+}=A^j_{i+}\frac {w^j_{i}-2u_2}{v^j_{i+1}+2u_1}\frac {v^j_i+2u_1}{w^j_i+2u_2}\cdot$$ (47)

Figure 3: The analytic solution for u1/u2=4 and v0=0.25u1 for the case $F_\parallel =0$. The downstream and upstream distributions are shown by solid and dashed lines respectively. The height of the vertical lines corresponds to the coefficients in Eqs. (34) and (35).

It is easy to show that for large w^j_i and v^j_i the change of velocity

$$% w^j_i-v^j_i\approx v^j_{i+1}-w^j_i\approx \frac 23(u_1-u_2).$$ (48)

Using Eq. (47) we find that for high velocities

$$% A^j_{i+}/(v^j_i)^2\sim (v^j_i)^{-3u_1/(u_1-u_2)}, \ v^j_i \gg u_1,u_2$$ (49)

that is exactly the spectrum of particles produced by diffusive shock acceleration.

The solution for u₁/u₂=4 and v₀=0.25u₁ is shown in Fig. 3. The critical velocities are $v_{\rm c-}\approx 0.2u_1$ and $v_{\rm c+}\approx 0.35u_1$ in this case. The far upstream distribution is shown by the dashed line on the left. The forward moving particles of this distribution reach the downstream region with slightly higher velocity (solid line on the left). The backward moving particles are directly accelerated along the characteristic and return upstream with a velocity close to 8u₁/3 (first dashed line on the right). Particles are accelerated further, moving between downstream and upstream regions of the shock (the other lines on the right). The similar $\delta$-function spectrum was obtained by Bogdan & Webb (1987) for cosmic ray acceleration in the so-called two-streaming approximation of Fisk & Axford (1969).

The spectrum will be smoother if one takes into account the scattering of particles in the transition region.

5 Numeric modeling of collisionless shocks

We use the system of Eqs. (12), (19), (22), (23) for the modeling of steady-state one-dimensional nonrelativistic shocks. We neglect the dynamical effects of the magnetic field, which means that the magnetic field produces only kinematic effects that are essential for oblique shocks. The upstream plasma state is then determined only by the sonic shock Mach number M and by the ratio of electron and ion temperatures  $T_{\rm e}/T_{\rm i}$.

A so-called Alfvén heating is not included in our model. It is well known that the energetic particles upstream of the shock can effectively generate Alfvén waves due to streaming instability (Wentzel 1969; Bell 1978). The damping of these waves results in the heating of the shock precursor (Völk & McKenzie 1982; McKenzie & Völk 1982). This heating is essential, since it decreases the effective Mach number and compression ratio of the shock. As shown in the Appendix, this effect can be formally included using the relation between sonic and Alfvén Mach numbers of the shock:

$$% M^2=\frac {M_{\rm s}^2M_{\rm a}}{M_{\rm a}+\frac {20}{39}M_{\rm s}^2}$$ (50)

The results obtained below for shocks with sonic Mach number M without Alfvén heating can be formally used for shocks with Alfvén heating and with Alfvén and sonic Mach numbers $M_{\rm a}$ and $M_{\rm s}$ respectively.

The details of the numeric method are given in the Appendix. The plasma flow with unity velocity and unity density enters the simulation box from its left boundary. The pressure at the right boundary was adjusted according to the motion of the shock in order to reach the steady state. We used the free escape boundary condition for the particle distribution at the left boundary.

We prescribe the dependence of the free path of particles on the velocity v' and space coordinate z. We used the free path $\lambda$ that is independent on z. This is not a strong limitation. Indeed, in the rather general case of separable dependence of $\lambda$on v' and z, the results obtained depend only on the integral $\int {\rm d}z/\lambda (z)$. Thus our results can be used also for this more general case.

The energy dependence of the mean free path $\lambda$ is determined by the spectrum of magnetic inhomogeneities and by the mechanism of scattering in the vicinity of a pitch angle of 90 degrees. Since the low frequency waves propagating along the magnetic field are subject to nonlinear wave steepening, one can expect that the magnetic spectra are close to k^-2 and the corresponding free path does not depend on energy. Such spectra (or slightly flatter ones) were indeed observed in the hybrid plasma simulations of parallel collisionless shocks (Giacalone et al. 1993; Giacalone 2004) and in the vicinity of the Earth bow shock in the solar wind. In the results presented here we used the energy independent free path.

We limit ourselves to the case of cold electrons $T_{\rm e}/T_{\rm i}=0$, as the dissipation of the Alfvén waves in the precursor may result in the preferable heating of ions. The solar corona is an example of this.

We start with the case of parallel shocks with Mach numbers 3.87 and 7.75. According to Eq. (50) the corresponding Alfvén Mach numbers are 7.5 and 30 that can be applied to the interplanetary and supernova shocks respectively. The downstream particle distributions are shown in the top panels of Figs. 4 and 5. The velocity and gas pressure profiles are shown in the bottom panels.

Figure 4: The results of simulation of the parallel shock with Mach number 3.77, $T_{\rm e}/T_{\rm i}=0$ and energy independent free path $\lambda =0.2$. The forward and backward particle distributions downstream N+ and N- are shown on the top panel by solid and dashed lines respectively. Plasma velocity u, pressure of ions P and electric force $F=qE_\parallel$are shown on the bottom panel by solid, dashed and dotted lines respectively. The total compression ratio r=4.05 is obtained.

Figure 5: The results of simulation of the parallel shock with Mach number 7.75, $T_{\rm e}/T_{\rm i}=0$ and energy independent free path $\lambda =0.2$. The forward and backward particle distributions downstream N+ and N- are shown on the top panel by solid and dashed lines respectively. Plasma velocity u, pressure of ions P and electric force $F=qE_\parallel$are shown on the bottom panel by solid, dashed and dotted lines respectively. The total compression ratio r=6.67 is obtained.

The simulated shocks demonstrate clear sub-shocks with a width of about one half of the free path $\lambda$ and a smooth precursor produced by the particles in the high energy tail that is seen in the top panels of Figs. 4 and 5. It was found that the main part of the downstream distribution in Figs. 4 and 5 is the adiabatically compressed far upstream distribution which was taken to be Maxwellian. The total electric potential corresponding to the electric force  $qE_\parallel$ is about 0.2-0.3 mu₁² according to Figs. 4 and 5.

The total compression ratio of the shock can be found from the following simple estimate.

In the precursor region the thermal particles are heated adiabatically. When their velocities become comparable to the plasma velocity, the backward moving particles are accelerated more efficiently (see Fig. 2) and injected for further acceleration. The compression ratio should be large enough for this to be the case. For estimation we assume that distribution of upstream particles is proportional to  $\delta (v'-v_{{\rm T}})$, where  $v_{\rm T}=\sqrt{3/5}u_1/M$ is the thermal velocity of the plasma. The maximum shock downstream velocity u₂ that is necessary for injection of such particles according to Fig. 2 and characteristics (26) is

$$% u_2=\frac 12v_{\rm T}\left( 8\frac {u_1}{v_{\rm T}}-3\right) ^{1/4}.$$ (51)

Thus the thermal velocity $v_{\rm T}$ should be larger than the critical velocity  $v_{\rm c-}$ (see Sect. 3). The second term in parentheses can be neglected. Writing the thermal velocity $v_{\rm T}$ in terms of Mach number M we obtain:

$$% r=2^{1/4}\left( \sqrt{5/3}M\right) ^{3/4}\approx 1.44M^{3/4}.$$ (52)

This formula is in good agreement with the simulated compression ratio. It was found that it is valid also for smaller free paths, when the maximum energy of the high energy tail is larger.

Figure 6: Comparison of our model with the cosmic ray modified shock approach for a parallel shock with Mach number 7.75, $T_{\rm e}/T_{\rm i}=0$ and energy independent free path $\lambda =0.2$. The downstream distributions for our model and for the cosmic ray approach are shown in the top panel by solid and dashed lines respectively. Plasma velocity u from our simulations, pressure of ions P and velocity profile for the cosmic ray approach are shown in the bottom panel by dotted, dashed and solid lines respectively.

This compression ratio is close to the maximum possible compression ratio of the idealized system that includes the infinitely thin gas sub-shock and the viscous precursor produced by energetic particles. Using the Ranke-Hugoniot conditions at the thermal sub-shock one can obtain for this case:

$$% r=2.5\left( M^2/5 \right) ^{3/8}\approx 1.37M^{3/4}.$$ (53)

This number is also close to the value found in numerical simulations of shock waves, modified by the cosmic ray pressure (see e.g. Berezhko & Ellison 1999). The compression ratio of the thermal sub-shock is 2.5.

The comparison of our results and results obtained for this simplified description are shown in Fig. 6. The Eq. (14) was solved in the upstream region together with the Euler equations. The top panel shows the comparison of downstream spectra obtained in these approaches. The velocity profiles are compared in the bottom panel. We find a very good agreement of these two methods. We estimate the injection efficiency at injection velocity $v_{{\rm inj}}=2u_{{\rm sub}}$to be about 0.024. Here $u_{{\rm sub}}$ is the sub-shock velocity. This number is in reasonable agreement with results of hybrid simulations (cf. Giacalone et al. 1993; Ellison et al. 1993) and the Earth bow shock observations (cf. Ellison et al. 1990).

We have also performed modeling of a quasi-parallel shock. The results obtained for a Mach number of 7.75 and an angle between the shock normal and magnetic field of $\theta =15^{\circ }$ are shown in Fig. 7.

Figure 7: Simulation of the quasi-parallel shock with $\theta =15^{\circ }$, Mach number M= 7.75, $T_{\rm e}/T_{\rm i}=0$ and energy independent free path $\lambda =0.2$. The forward and backward particle distributions N+ and N- are shown in the top panel by solid and dashed lines respectively. Plasma velocity u, pressure of ions P and electric force $F=qE_\parallel$shown in the bottom panel by solid, dashed and dotted lines respectively. The total compression ratio r=7.21 is obtained.

For such oblique shocks the dip appears between the main part of the thermal distribution and the high-energy tail (see the top panel). The injection efficiency is the same as in the previous case. $15^{\circ}$ is the maximal value of $\theta$ that allows us to model steady state shocks with M=7.75. Even this value is rather large, because it corresponds to the normal angle 40$^{\circ}$ upstream of the thermal subshock. For larger $\theta$ this dip becomes more pronounced and the thermal subshock width drops to a grid step. In some cases some instability in the downstream region was observed. It is possible to model the larger values of $\theta$ for the smaller Mach numbers (e.g. $\theta =30^{\circ}$ for M=3.87). Thus our approach can be used only for quasiparallel shocks with the normal angle less than 40$^{\circ}$ just upstream of the thermal subshock.

If the magnetic field is oblique enough, it prevents particles from returning from downstream along magnetic lines and the heating of particles is impossible. Since it is necessary to have the energetic particles downstream for the existence of a shock, we expect that the transition region will be squeezed to widths of the order of the ion gyroradius. Heating then may take place in a different regime, where the drift motion of particles perpendicular to magnetic field is essential. The injection efficiency also may be very different.

6 Conclusion

Many astrophysical problems deal with collisionless plasma. Different low frequency instabilities may produce magnetic inhomogeneities that can scatter particles and play the role of thermal collisions. When unstable waves propagate preferably along the mean magnetic field, the resonant scattering of particles is weak in the vicinity of 90 degrees pitch angle. In this case it is possible that particle distributions are almost isotropic in both hemispheres.

Following the approach of Isenberg (1997) we derive the general transport Eq. (12) (or equivalent Eq. (9)) which takes into account the motion of the medium and energy changes of particles. Our consideration is not limited by high particle velocities and by small anisotropies of particle distribution.

In the case of high particle velocities and small anisotropies the equation derived can be reduced to the standard cosmic ray transport equation. In the case of negligible energy changes and advective transport the equation may be transformed to the so-called telegraph equation.

The equations derived can be solved together with magnetohydrodynamic equations (see Sect. 3). We used the fluid approximation for the electron transport. This is justified if electrons are effectively scattered by magnetic inhomogeneities. If this is not the case, one can use Eqs. (12) for electron transport. In this way electron injection at collisionless shocks can be investigated.

The transport Eqs. (12) with Eqs. (19), (22) and (23) can be used for the solution of different astrophysical problems when the self-consistent determination of the flow velocity, magnetic field etc. is necessary.

We apply the equation derived to investigate particle acceleration and injections at astrophysical shocks. For the prescribed velocity profile we have found the exact analytical solution (see Sect. 4).

We also have performed the modeling of collisionless shocks and found a formula for the shock compression ratio (see Eq. (52)). This ratio is very close to the maximal possible value in the framework of the simplified approach used for cosmic ray modified shocks. The corresponding thermal sub-shock ratio is 2.5.

We found that this feature is universal and does not depend on the maximum energy of the accelerated particles. This was checked in our simulations up to the maximum momentum $p_{\max }\sim 100 \ mu_1$ which was limited by the step of our uniform spatial grid. This property will not change for higher maximum energies corresponding to the supernova environment. On the other hand we found that the injection efficiency $\eta$ is not universal. For $p_{\max }\sim 10\ mu_1$ it is about 0.024 for an injection velocity equal to two sub-shock velocities. It slowly decreased as the maximum energy increases and is about 0.003 for young supernova remnants. This is because at small maximum energies the thermal ions "feel'' a larger compression ratio than the thermal subshock compression ratio since the subshock and precursor width are not strongly different.

We found that the equations derived can be used to investigate injection and acceleration in quasiparallel shocks with the subshock normal angle $\theta <40^{\circ}$. The dissipation at more oblique shocks is described by a different mechanism which takes into account the motion of particles in the direction perpendicular to the magnetic field. The injection efficiency at such shocks is unknown, but it is possible that it is rather low. We confirm the previous findings that quasiparallel shocks are very efficient accelerators with high injection efficiency (see e.g. Ellison & Eichler 1984; Giacalone et al. 1993; Malkov & Drury 2001) which does not depend on the angle $\theta$ (Giacalone et al. 1997).

Acknowledgements

The author is grateful to Heinz Völk for fruitful discussions of the collisionless shock physics. The author also thanks the anonymous referee for important comments.

References

 
• Achterberg, A. 1981, A&A, 98, 161  
• Bell, A. 1978, MNRAS, 182, 147  
• Berezhko, E. G., & Ellison, D. C. 1997, ApJ, 526, 385  
• Bogdan, T. J., & Webb, G. M. 1987, MNRAS, 229, 41  
• Dolginov, A. Z., & Toptygin, I. N. 1966, JETP, 51, 1771  
• Gleeson, L. J., & Axford, W. I. 1969, ApJ, 149, L115  
• Earl, J. A. 1974, ApJ, 188, 379  
• Ellison, D. C., & Eichler, D. 1984, ApJ, 286, 691  
• Ellison, D. C., Möbius, E., & Paschmann, G. 1990, ApJ, 352, 376  
• Ellison, D. C., Giacalone, J., Burgess, D., & Schwartz, S. J. 1993, JGR, 98, 21085  
• Fisk, L. A., & Axford, W. I. 1969, Sol. Phys., 7, 486  
• Fisk, L. A., Schwadron, N. A., & Gloeckler, G. 1997, Geoph. Res. Lett., 24, 93  
• Giacalone, J. 2004, ApJ, 609, 452  
• Giacalone, J., Burgess, D., Schwartz, S. J., & Ellison, D. C. 1993, ApJ, 402, 550  
• Giacalone, J., Burgess, D., Schwartz, S. J., Ellison, D. C., & Benett, L. 1997, JGR, 102, 19789  
• Godunov, S. K. 1971, in Equations of Mathematical Physics, Nauka, Moscow (in Russian)  
• Gombosi, T. I., Jokipii, J. R., Kota, J., & Williams, L. L. 1993, ApJ, 403, 377  
• Isenberg, P. 1997, JGR, 102, 4719  
• Jokipii, R. 1966, ApJ, 146, 480  
• Kota, J. 2000, JGR, 105, 2403  
• Kota, J., Jokipii, J. R., Kopriva, D. A., et al. 1982, ApJ, 254, 398  
• Krymsky, G. F. 1964, Geomagn. Aéron., 4, 977  
• Kulsrud, R. M. 1983, in Basic Plasma Physics I, ed. A. A. Galeev, & R. N. Sudan, North-Holland, New York, 115  
• Leroy, M. M., & Winske, D. 1983, Ann. Geophys., 1, 527  
• Malkov, M. A., & Drury, L. O'C. 2001, Rep. Progr. Phys., 64, 429  
• McKenzie, J. F., & Völk, H. J. 1982, A&A, 116, 191  
• Quest, K. B. 1986, JGR, 91, 8805  
• Quest, K. B. 1988, JGR, 93, 11311  
• Parker, E. N. 1961, J. Nucl. Energy, C2, 146  
• Parker, E. N. 1965, Planet. Space Sci., 13, 9  
• Skilling, J. 1971, ApJ, 170, 265  
• Völk, H. J. 1973, Astroph. Space Sci., 25, 471  
• Völk, H. J. 1975, Rev. Geophys. Space Phys., 13, 547  
• Völk, H. J., & McKenzie, J. F. 1982, Proc. 17th ICRC, Paris, 9, 246  
• Wentzel, D. G. 1969, ApJ, 156, 303

 

  
7 Online Material

Appendix A: Formal inclusion of Alfvén heating

The heating of the cosmic ray precursor is described by the following equations (Völk & McKenzie 1982; McKenzie & Völk 1982):

$$% \frac {\partial }{\partial z}u\rho =0$$ (A.1)

$$% \rho u\frac {\partial u}{\partial z}=-\frac {\partial P_{\rm g}}{\partial z} -\frac {\partial P_{\rm c}}{\partial z},$$ (A.2)

$$% u\frac {\partial P_{\rm g}}{\partial z}+\gamma _{\rm g}P_{\... ...{\rm g}-1)V_{\rm a}\frac {\partial P_{\rm c}}{\partial z}\cdot$$ (A.3)

Here $V_{\rm a}$ is the component of Alfvén velocity parallel to the shock normal, $P_{\rm g}$ and $P_{\rm c}$ are the pressure of the gas and cosmic rays respectively, and  $\gamma _{\rm g}$ is the gas adiabatic index. The third equation describes the gas heating due to the damping of Alfvén waves generated by the cosmic ray streaming instability. The second equation is the Euler equation of motion.

The gradient of the cosmic ray pressure can be found from Eq. (A.2) and substituted into Eq. (A.3). For high Mach number shocks the solution can be written as

$$% P_{\rm g}=\left( P_{\rm g1}+\rho _1u_1V_{\rm a1}\frac {\gamma _... ... _{\rm g}-1}{\gamma _{\rm g}+\frac 12} \left( \frac {u}{u_1}\right) ^{1/2}\cdot$$ (A.4)

Here $P_{\rm g1}$ and $V_{\rm a1}$ are the gas pressure and Alfvén velocity in the medium, in which the shock propagates with speed u₁.

For high Mach number shocks the Alfvén heating is essential only at the very beginning of the precursor. In the rest, the gas is heated adiabatically and the last term in Eq. (A.4) can be neglected. The gas pressure  $P_{\rm g1}$and Alfvén velocity  $V_{\rm a1}$ can be expressed in terms of Alfvén and sonic Mach numbers $M_{\rm a}$ and $M_{\rm s}$ respectively. As a result, we obtain the sonic Mach number M of the shock without Alfvén heating, which is equivalent to the shock with the Alfvén heating:

$$% M^2=\frac {M_{\rm s}^2M_{\rm a}}{M_{\rm a}+\gamma _{\rm g}\... ...{\gamma _{\rm g}-1}{\gamma _{\rm g}+\frac 12}M_{\rm s}^2}\cdot$$ (A.5)

For $\gamma _{\rm g}=5/3$ we obtain the formula (50) from the main part.

Appendix B: Numeric method

We used the implicit finite difference scheme for the solution of the one-dimensional version of Eq. (12):

$$% \frac {\partial n}{\partial t}=L_z n+L_vn.$$ (B.1)

Here n describes the pair of forward and backward particle distributions $n_\pm =v'^3N_\pm$, and L_z and L_v are the operators in the coordinate and velocity space respectively:

$$% L_zn=-\frac {\partial }{\partial z}\left( u\pm \frac 12v'b_z\right) n_\pm,$$ (B.2)

$$% L_vn=-v'\frac {\partial }{\partial v'}\left( \pm \frac {F_\... ... u}{\partial z}\right) n_\pm + \nu _\mp n_\mp -\nu _\pm n_\pm.$$ (B.3)

The finite difference version of Eq. (B.1) was solved in two steps (Godunov 1971):

$$% \frac {n_{k+1/2}^{i,j}-n_k^{i,j}}{\tau }=L^{\rm dif}_vn_{k+1/2}^{i,j}+L^{\rm dif}_{z}n_{k}^{i,j},$$ (B.4)

$$% \frac {n_{k+1}^{i,j}-n_{k+1/2}^{i,j}}{\tau }=L^{\rm dif}_zn_{k+1}^{i,j}-L^{\rm dif}_{z}n_{k}^{i,j}.$$ (B.5)

Here n_k^i,j is the value of distribution n(z,v',t) at $z=i\Delta z$, v'=v_j, $t=k\tau$, where $\tau$ is the time step, $\Delta z$ is the grid size in z-direction and v_j is the value of velocity v' at the knot of the velocity grid with number j. We use the following finite difference analog  $L^{\rm dif}_z$ of the operator $L_zn=-\partial /\partial zAn$ (Kota et al. 1982):

$$% L^{\rm dif}_zn_k^{i,j} \left( \frac 16A_k^{i+2}n_k^{i+2,j}-A_k^{i+1}n_k^{i+1,j}+ \frac 12A_k^{i}n_k^{i,j} \right.$$ =

$$\left. + \frac 13A_k^{i-1}n_k^{i-1,j}\right) /{\Delta z}, \ A_k^i<0,$$ (B.6)

$$% L^{\rm dif}_zn_k^{i,j} \left( -\frac 16A_k^{i-2}n_k^{i-2,j}+A_k^{i-1,j}n_k^{i-1,j}- \frac 12A_k^{i}n_k^{i,j} \right.$$ =

$$\left. - \frac 13A_k^{i+1}n_k^{i+1,j}\right) /{\Delta z}, \ A_k^i>0.$$ (B.7)

A similar operator was used for the finite difference analog  $L_v^{\rm dif}$ of the velocity operator (B.3).

The numeric scheme (B.4), (B.5) with these operators has a third order accuracy on z and v' and first order accuracy on t.

The parallel force $F_\parallel$ was recalculated at each time step according to Eq. (21). The medium velocity u_kⁱ at each time step was calculated according to the finite difference version of the Euler Eq. (19):

$$% \frac {\rho _{k+1}^iu_{k+1}^i-\rho _{k}^iu_{k}^i}{\tau } -\frac {\rho _{k+1}^{i}(u_{k+1}^{i})^2-\rho _{k+1}^{i-1}(u_{k+1}^{i-1})^2}{\Delta z}$$ =

$$+ \frac {P_{k+1}^{i-1}-P_{k+1}^{i+1}}{2\Delta z}.$$ (B.8)

Here P_kⁱ is the sum of the ion and electron pressures. The ion pressure and plasma density were calculated from the particle distribution at the each time step. The electron pressure was found from the adiabatic equation of state. Equation (B.8) is the quadratic equation for the velocity u_k+1ⁱ. The gas pressure at the right boundary  $P_k^{\rm b}$ was prescribed in accordance with the shock motion:

$$% P_k^{\rm b}=P_{k-1}^{\rm b}+\frac 12(Z_k^{\rm s}-Z_{k-1}^{\rm s}).$$ (B.9)

Here $Z_k^{\rm s}$ is the position of the shock. When the shock moves to the right boundary, the pressure increases. This permits us to reach steady state.

We use the hyperbolic tangent for the initial velocity profile. The initial particle distribution was Maxwellian with a temperature and density dependence corresponding to the momentum and mass conservation. At times of about 60-100 in dimensionless units the system reaches the steady state. We use 200 grid points in the z-direction, 100 grid points for the logarithmic v' grid and the time step between 0.005 and 0.1 for different runs. We obtained the total momentum, energy and mass conservation with several percent accuracy.