A kinetic model for runaway electrons in the ionosphere

Abstract. Electrodynamic models and measurements with satellites and incoherent scatter radars predict large field aligned current densities on one side of the auroral arcs. Different authors and different kinds of studies (experimental or modeling) agree that the current density can reach up to hundreds of µA/m 2 . This large current density could be the cause of many phenomena such as tall red rays or triggering of unstable ion acoustic waves. In the present paper, we consider the issue of electrons moving through an ionospheric gas of positive ions and neutrals under the influence of a static electric field. We develop a kinetic model of collisions including electrons/electrons, electrons/ions and electrons/neutrals collisions. We use a Fokker-Planck approach to describe binary collisions between charged particles with a long-range interaction. We present the essential elements of this collision operator: the Langevin equation for electrons/ions and electrons/electrons collisions and the Monte-Carlo and null collision methods for electrons/neutrals collisions. A computational example is given illustrating the approach to equilibrium and the impact of the different terms (electrons/electrons and electrons/ions collisions on the one hand and electrons/neutrals collisions on the other hand). Then, a parallel electric field is applied in a new sample run. In this run, the electrons move in the z direction parallel to the electric field. The first results show that all the electron distribution functions are non-Maxwellian. Furthermore, runaway electrons can carry a significant part of the total current density, up to 20% of the total current density.


Introduction
The existence of large field-aligned current densities in narrow auroral structures has been inferred over the last years by Correspondence to: G. Garcia (gga@cetp.ipsl.fr)using satellites (Cerisier et al., 1987;Berthelier et al., 1988), incoherent scatter radars (Rietveld et al., 1991) and numerical models (Noël et al., 2000).Cerisier et al. (1987) interpreted a magnetic pulse recorded by the magnetometer on board the AUREOL 3 low altitude satellite as the signature of current densities as high as 500 µ A/m 2 in a current sheet 20 m wide.Later on, Stasiewicz et al. (1996) reported smallscale current densities observed on board the Freja satellite of a few hundred µ A/m 2 .Recently the measurements made from the ØRSTED satellite has enabled the detection of finescale structures, as low as 75 m, in the high latitude fieldaligned current system.Very intense but thin sheets or narrow filaments of field-aligned currents (FAC) up to several hundreds of µ A/m 2 have been reported by Stauning et al. (2003).At higher altitudes of 1000 to 4000 km, the downward Birkeland currents are carried by suprathermal electrons at energies from 10 to 500 eV and fluxes greater than 10 9 electrons.cm−2 .s−1 (Klumpar and Heikkila, 1982;Carlson et al., 1998).Klumpar and Heikkila (1982) suggested that they are runaway electrons from the ionosphere produced by a downward field-aligned component of the electric field.
Radar observations have also suggested the existence of extremely intense current densities.A large increase in the electron temperature measured in filamentary aurora with the European incoherent scatter radar has also been interpreted as a hint of the presence of intense FAC densities (Lanchester et al., 2001).They modeled the observations with an 1-D electron transport and ion chemistry code.They concluded that, to account for the observed changes in the electron temperature, a source of electron heating was required in addition to local heating from energy degradation of precipitating electrons.They showed that Joule heating in a strong FAC of 400 µ A/m 2 can account for the required heat source.Strong enhancement in incoherent scatter radar spectra have also been observed by Rietveld et al. (1991) and interpreted as unstable ion-acoustic waves triggered by large FAC densities.Threshold calculations for the two-stream instability for typical ionospheric parameters lead to FAC densities carried G. Garcia and F. Forme: A kinetic model for runaway electrons by thermal electrons which have to be in excess of 1000 µ A/m 2 (Rietveld et al., 1991).St Maurice et al. (1996) suggested that large parallel current densities carried by thermal electrons can be triggered in the ionosphere with horizontal scale sizes of a hundred meters or less.
In order to explain the large FAC in the ionosphere, Noël et al. (2000) has developed a two-dimensional model of short-scale auroral electrodynamics that uses current continuity, Ohms's law, and 8-moment transport equations for the ions and electrons in the presence of large ambient electric field to describe wide auroral arcs with sharp edges in response to sharp cut-offs in precipitation.Using an ambient perpendicular electric field of 100 mV/m away from the arc and for which electron precipitation cuts off over a region 100 m wide, they showed that parallel current densities of several hundred µ A/m 2 can be triggered together with a parallel electric field of the order of 0.1 mV/m around 130 km altitude.In a rather similar model, Otto et al. (2003) showed that ohmic heating due to intense FAC densities up to 600 µ A/m 2 can lead to the formation of tall red rays.The resulting heating leads to an electron temperature in excess of 10000 K in the upper F-region.
Therefore if one wants to understand the electrodynamics of the auroral arcs, the role of the ionosphere in the generation of intense parallel currents and the associated parallel electric fields is a matter of particular interest.However, it is well-known that if an electric field (not too weak) is applied to a collisional plasma, some electrons experience unlimited "runaway" acceleration (Dreicer, 1959).The reason is straightforward and well-known: the friction force acting on an electron travelling with velocity v is a non-monotonic function, having a global maximum at the thermal speed.For an electron moving faster than this speed, the collision frequency decreases with increasing velocity.Therefore, if a sufficiently fast electron starts accelerating in the electric field, the dynamical friction force decreases.A critical electric field, known as the Dreicer field E c , has been calculated by Dreicer (1959).It is a measure of the electric field which is required if the drift velocities are to increase and exceed the most probable random speed in one free collision time.However, using the Dreicer field is only an estimate of the importance of kinetic effects because runaway still occurs for E<E c (Dreicer, 1959).
The acceleration of runaway electrons have been studied in: solar flares (e.g.Moghaddam-Taaheri andGoertz, 1990), Tokamaks (e.g. Liu et al., 1977) and red sprites (e.g.Bell et al., 1995).In the ionosphere, the large field-aligned current densities have only been modeled using fluid models (Noël et al., 2000;Otto et al., 2003;Noël et al., 2005).However, the fluid models could be altered by runaway electrons.Using typical ionospheric parameters, Otto et al. (2003) estimated the Dreicer field to be of about E c ≈4×10 −5 V/m which is much higher than typical F-region electric fields, from Ohm's law, of about 10 −6 V/m.On the other hand, Noël et al. (2000) published a parallel electric field of about 5×10 −4 V/m in the E region suggesting that a substantial part of the electron distribution function could be freely accelerated.Papadopoulos (1977) suggested that runaway electrons from intense FAC could create non-Maxwellian electron distribution functions (EDF) that could in turn trigger Langmuir turbulence in the ionosphere.In their hypothesis, the thermal ionospheric electrons were accelerated by a parallel electric field due to anomalous resistivity.However, a quantitative estimate of runaway electrons in the low altitude ionosphere has never been done.
This study is a first step of a kinetic model of the highly collisional low altitude ionosphere (E and low F-region).In Sect.2, we will mainly focus on the description and the tests of the collision operators which is of a crucial importance for this study .Then, in Sect.3, we will consider the issue of electrons moving through an ionospheric gas of positive ions and neutrals under the influence of a static electric field similar to Noël et al. (2000).

Description of the kinetic model
First, we are interested in charged particle collisions.In this section we will present the methods that are used to simulate these collisions.We want to investigate the interactions between charged particles in a highly collisional plasma.For this purpose, we use a Fokker-Planck approach (see Sect. 2.1.1),which describes binary collisions between charged particles with long-range interactions.Effectively, in our case, the long-range interactions are more dominant than the short-range interactions, since the coulomb logarithm ln = ln(λ d /p 0 ) 1, where λ d is the Debye length and p 0 is the impact parameter (Rosenbluth et al., 1957).Under typical ionospheric conditions, ln is around 15.We also consider that the collisions are binary as the impact parameter (p 0 ≈4×10 −19 m) is smaller than the mean length between particles (d e ≈2×10 −10 m).

The Fokker-Planck approach
We consider charged particles of species a interacting with species b.Species a are electrons and species b are either electrons or ions.All the particles collide, but the collisions with long-range interactions, which correspond to small pitch angle, play a more dominant role than those at a closer range.
The classical Fokker-Planck equation including electron/electron (e/e) and electron/ion (e/i) collisions (Krall and Trivelpiece, 1986) is: friction diffusion where f a , the electron distribution function (EDF), is a function of position, velocity and time, v a is the electron velocity and F ext are the external forces (in our case the electric force).
The friction and diffusion coefficients on the right hand side of equation 1 are used in the Langevin equation (see Sect. 2.1.2).These coefficients have been given by Rosenbluth et al. (1957) and MacDonald et al. (1957): where: The coefficients G ab and H ab , which govern the diffusion and dynamic friction, are scalar functions of the vector velocity v: where: m a is the electron mass, m b is the mass of the target particles, e is the elementary charge, v a is the electron velocity, v b is the velocity of target particles and f (v b ) is the velocity distribution function of the target particles.
If the distribution functions of target particles are isotropic, only three coefficients need to be considered (Manheimer et al., 1997): the friction coefficient The assumption of isotropic scatterers is a reasonably good approximation, since all electron collisions tend to isotropize the EDF.The thermal part of the EDF isotropizes particularly rapidly, and e-e scattering of any electron (even a fast one) is normally dominated by scattering off thermal electrons.Manheimer et al. (1997) considered the question of the limits of validity of this approximation for very anisotropic situations.They conclude that this approximation retains quantitative accuracy in situations where the EDF is single-peaked.We can infer that: where: D 11 is the half angular diffusion coefficient and D 33 is the longitudinal diffusion coefficient The analytical expressions for the diffusion coefficients with Maxwellian distribution functions are well known (Barghouthi and Barakat, 2005 and references therein).However, since we expect a non-Maxwellian EDF, we perform a numerical integration of these coefficients at each time step.
In the following, we also assume that: -The ion distribution function is a stationary non-drifting Maxwellian, since the i/e relaxation time is longer than the e/i relaxation time (τ i/e ≈10 4 τ e/ i ).
-The plasma is considered to be quasi-neutral.

The Langevin equation
In order to use the collision operators defined in the previous section and to calculate the new velocity for each electron at each time step, it is necessary to go from the Fokker-Planck equation to the Langevin equation.This equation is equivalent to the Fokker-Planck equation.To first order accuracy in t, the 3-D Langevin equation takes the following form: where: is the friction coefficient, F ext is the electric force, m is the electron mass and Q= is a random velocity vector that corresponds to the variation of the velocity in www.ann-geophys.net/24/2391/2006/Ann.Geophys., 24, 2391-2401, 2006 G. Garcia and F. Forme: A kinetic model for runaway electrons the three direction due to the diffusion.It is chosen from the distribution (Manheimer et al., 1997): where D 11 corresponds to 1 2 (Eq. 9) and D 33 to (Eq. 10).The friction tends to slow down the electrons whereas Q scatters them.

The conservation of energy and linear momentum
The Fokker-Planck equation ( 1) conserves momentum and energy.However, the numerical implementation of e/e and e/i scattering can lead to small deviation from energy and momentum conservation.The diffusion coefficient of the Langevin equation implies the choice of a random velocity Q in a distribution.On average, these increments will conserve the energy and momentum.However, we use a finite number of particles N and as a consequence an error of √ N can always occur.Since the drift velocity is small compared to the thermal velocity, the error is negligible.However, we make corrections to restore the conservation laws.We calculate the energy before and after the scattering, and then, we renormalize the velocity for each electron (Manheimer et al., 1997): where W j is the total kinetic energy of the electrons before the collisions and W j is the the total kinetic energy after the collisions.

The electron/neutral collisions: Monte-Carlo method
Since we want to model collisions in the ionospheric Eregion, we must include electron/neutral (e/n) collisions.The simulation procedure must give us the time interval between each pair of collision and the change in the electron velocity due to each collisions (Winkler et al., 1992).For a collision with a collision frequency ν that is independent of the relative velocity, the probability that a test particule will suffer no collision for a time interval t is given by: If the collision frequencies were constant, the time interval between collisions would be generated by using: where R is a random number with a uniform probability between 0 and 1 (Lin and Bardsley, 1977).However for the collisions that we are considering, the total cross section of interaction is speed-dependent, and the collision frequency is consequently a function of energy.The time interval therefore depends on the continually changing relative velocity of the colliding particles, so that the simple expression for the probability that the time between two collisions has the value t given by Eq. ( 15) is no longer valid.To correct the speed dependence difficulty, we use a "null collision" approach, described in the next section and a Monte-Carlo method.

Time of free flight: "null collision" approach
The time of free flight t f (Skullerud, 1968;Lin and Bardsley, 1977) is: where: r f is a random number chosen from a uniform distribution in the range between 0 and 1 (i.e.[0,1]) -and, where, ν is the total e/n collision frequency (elastic and inelastic) and ν null is the null collision frequency.The null collision frequency is chosen in such a way to keep ν tot constant, i.e. ν tot is the maximal collision frequency in the considered energy range.
The classical ν tot is the sum of all the collision frequencies e.g.elastic and inelastic collision frequencies.We introduce a null collision frequency since the electron velocity varies as a function of time.

Types of collisions
The probability of each collision is: with P coll el + P coll in + P coll null = 1, where ν el , ν in and ν null are the elastic, inelastic and total collision frequency, respectively.
The P coll el , the P coll in and the P coll null are the probabilities to have an elastic, inelastic and no collision, respectively.
We model the elastic and inelastic collision frequency: where n n is the neutral density, σ el , σ in are the tabulated elastic and inelastic collision cross-sections and v e is the electron velocity.
Elastic and inelastic cross-section as a function of electron velocity are given by Hubner et al. (1992).We assume that the neutrals are at rest.When an electron undergoes an elastic collision, its energy is conserved.So, we only diffuse the electron velocity in space (Yousfi et al., 1994): v y = v(sin χ sin η cos φ + sin χ cos η cos θ sin φ where: -η=2πr η where r η is a random number chosen between 0 and 1.This results in η being uniformly distributed in the range of [0, 2π ].
-The angle of deflection χ is assumed to be isotropic and chosen in the range [0, π] -The angles θ and φ are the polar and azimuthal angle in the laboratory frame of the incident velocity vector, v, respectively.
For an inelastic collision, the electrons will lose some of their energy.To determine the amount of lost energy, we tabulate the most likely reactions.Each reaction is energy dependent, therefore a knowledge of the energies that are involved is required.For loss of energy, we can calculate the corresponding loss of velocity, v l .The new velocity of the electron will be: where: v new is the electron velocity after the collision v old is the electron velocity before the collision v l is the velocity loss due to inelastic collision.It is tabulated in Gerjuoy and Stein (1955) and Gilmore (1965).
The electrons are then scattered in space using the same method as for the elastic collisions (23, 24, 25).

Test of the collision operator
In this section, we present numerical results as a test case of our collisions operators.The test involves the relaxation of an initially square distribution.We know that a square distribution has to relax towards a Maxwellian because of e/e and e/i collisions.We use n e =n i =5×10 10 m −3 , T e =T i =1500 K.We use 3×10 4 particles to represent the electrons.The electrons are scattered at intervals t=10 −4 ν −1 0 where ν 0 is the e/i collision frequency equal to 60 s −1 .

Time-evolution of the distribution function
In order to test the Fokker-Planck operator, the e/e and e/i collisions are the only physical processes considered in this simulation.We will study the time-evolution of an initial non-Maxellian rectangular EDF in velocity space.In the Fig. 1, we plot the frictional coefficients < v e > i < t> and < v e > e < t> in function of the electron velocity.The frictional term is much more important for weak velocities.The maximum of the friction coefficient is obtained when v=v thi =5.8×10 −3 v the .
In Fig. 2, panel a, we present the EDF as a function of time.At early times (0.1−0.5 ν −1 0 ), the distribution function becomes rounded but is still non-Maxwellian.Representing the logarithm of the distribution function as a function of the squared velocity, a Maxwellian distribution function would be represented by a straight line.The distribution does not change very much at low energy however the tail begins to sprawl at higher energy.The maximum of the squared velocity increases from 2.4 v 2 the over to 3 v 2 the between t=0.5 and t=10 ν −1 0 .This means that the maximum electron energy increases from 0.3 eV to 0.4 eV.The frictional and diffusion coefficients (Eqs.7 to 9) decrease rapidly with particle speed v (Fig. 1), so we expect the system to approach equilibrium faster for the electrons that are in the low-energy range, and slower for those in the high-energy tails.
At t=0.5 ν −1 0 , the EDF is close to the final Maxwellian in the thermal range but the high energy tail remains non-Maxwellian.At time, t=1 ν −1 0 , the EDF is Maxwellian over the entire energy range.In Fig. 2, panel b, we show that the total kinetic energy is conserved.This can be seen as a horizontal line as a function of time.
In Fig. 2, panel c, we present entropy as a function of time.We observe that the entropy increases strongly up to 0.6 ν −1 0 and then stabilizes.This is due to the fact that the modifications that are occurring in the tail of the distribution are too small to influence the entropy.
The second numerical simulation that we undertook was to examine the effects of the neutral collision operator on the EDF.This was achieved by only considering e/n collisions.We assumed the neutral species to be N 2 with a number density of 10 15 m −3 .The numerical results are shown in Fig. 3.One may see that the time to get a Maxwellian is longer than in the case of e/e and e/i collisions.We can consider that the relaxation time is 3 ν −1 0 , since the EDF is rounded and does not evolve significantly after t = 3 ν −1 0 (Fig. 3, panel a).The e/n collision frequency is ν −1 0 3 as a consequence, the relaxation time should be of the order of 3 ν −1 0 .We should also notice that the "final" distribution is isotropic but not Maxwellian.Indeed, the e/n collisions tend to isotropize the distribution but there are no frictional forces thereby preventing the EDF from becoming Maxwellian.
The energy panel (Fig. 3, panel b) shows that energy is no longer conserved.The inelastic collisions act to decrease the energy but the inelastic collision frequency is so small that the reduction of the total kinetic energy is weak.
The entropy (Fig. 3, panel c) also increases quickly up to 3 ν −1 0 .For times larger that this, the entropy increases slowly.So, this fact confirms that 3 ν −1 0 is the relaxation time to the equilibrium.
The results presented in this section show that our choice of collision operator is justified.Furthermore, the time scales seem very reasonable as the relaxation time for the e/i and e/n collisions are of the order of 1 ν −1 0 and a few ν −1 0 , respectively.We can also keep in mind that the loss of energy is due only to the inelastic e/n collisions.

Simulation model
Our goal was to investigate the dynamical behavior of the electron under the influence of an applied static electric field.At each time step, we calculate the EDF as a function of altitude.The parameters used in the model were chosen to be consistent with typical ionospheric values at 200 km altitude.Since our simulation is 1-D in space and 3-D in velocity, the lower boundary was chosen so that the perpendicular conductivities could be neglected.The perpendicular conductivities are dominant below 200-km altitude.We do not take into account the magnetic field since the electrons are strongly magnetized and can be regarded as firmly attached to a given magnetic field line.This is a good approximation as long as the electron Larmor radius is very small compared to any macroscopic length scale.We also neglect the perpendicular electric field that may exist at these altitudes.It is well known that when a perpendicular electric field, of more than 50 mV/m is present, the ion velocity distribution function (IDF) becomes non-Maxwellian (Hubert, 1982 andWinkler et al., 1992).As a first approximation, the non-Maxwellian IDF can be represented by a bi-Maxwellian.In order to quantify the possible effects on the coefficients D we have calculated the resulting friction and diffusion coefficients (not shown here).For velocities higher than the ion thermal velocity (about 10 −2 V the ) the coefficients remain the same.The major effect of a non-Maxwellian IDF is a decrease of the diffusion and the friction coefficients below 10 −2 V the .The friction coefficient being reduced, we can expect that more runaway electrons will be created.Therefore our results should underestimate the number of runaway electrons.A schematic representation of our model is presented in Fig. 4.
The 1-D simulation region of length (H =100 km) is divided into spatial cells of length h.The length of each cell was chosen so that it was less than a mean free path (h=2.5 km).The plasma in each cell is assumed to be uniform.Initially, particles are distributed randomly in each cell and in velocity space in terms of given plasma conditions such as density profile and temperature profile.We use T e =2000 K, T i =1000 K, the density is decreasing linearly as we move from the bottom of the box to the top.The electron and ion density varies from 10 11 to 5×10 10 m −3 while the neutral density varies from 2×10 15 to 3×10 13 m −3 .For the sake of simplicity we consider only one species of ion and neutral.The ion and neutral species being respectively O + and N 2 .We apply a parallel electric field (z aligned), which decreases from 3×10 −5 to 0 V/m .The electric field was chosen from Noël et al. (2000), their Fig. 6.
The electrons are scattered by electron/electron, electron/ion and electron/neutral collisions and the new velocity is calculated using the Langevin equation with a time step, t=10 −4 ν −1 0 (ν 0 =78 s −1 ).After the calculation of the new electron velocities, the electrons positions are found using: We only take into account motions along the z-direction.All the parameters are assumed to be homogenous in the perpendicular plane.We calculate the new EDF that needs to be taken into account in the Fokker-Planck equation.For each time step, the EDF is computed for each altitude and called f (z, v, t+ t).The EDF is obtained by computing the histogram for all velocities of the particles present in [z, z+h].We can then use the EDF to compute the friction and diffusion coefficients.In Fig. 5, we present a schematic of the computational scheme that was employed.

Boundary conditions
When a particle reaches the boundaries (z=0 or z=H) i.e. when it leaves the simulation region, we inject a new particle to conserve the overall density within the simulation box.We choose the velocity of the injected particle with the repartition function given by: In the case that the particle reaches z=H, the new particle is injected in the range [0, h].The new velocity is choosen according to G(v, z, t) where f (v, z, t) corresponds to f (v, 0, t).In the other case, where the particle reaches z=0, the new particle is injected in the range [H, H-h] with a new velocity chosen with G(v, z, t) where f (v, z, t) is f (v, H, t).The new altitude in the predetermined range is chosen randomly.

Results
In Fig. 6, we present the EDF for three different altitudes at different times (0, 1, 10, 25, 50 and 90 ν −1 0 ).The three different altitudes correspond to the bottom, middle and top of the simulation box (2.5, 50, 100 km).In Fig. 6, panel a, the drift velocity, defined as The EDF slowly shifts towards positive velocities as time increases.It may also be noticed that the slope of the EDF decreases as the electrons are heated by Ohmic disspation.The EDF remains symmetric with respect to the maximum of the EDF: each tail sprawls symmetrically.
As time increases, electrons from the low altitude subboxes can reach the middle altitudes.Fig. 6, panel b, shows the EDF at z=50 km.The drift velocities are v d =0 v the at t=0 ν −1 0 and v d =0.20 v the at t=90 ν −1 0 .This clearly shows that the electrons have been accelerated from the lower altitudes.In addition, the tails begin to sprawl as a function of time due to two effects: the bottom, middle and top 100 km).In 6, panel a, = vf (v)dv, increases from = 0.15 v the at t = 90 ν −1 0 .s positive velocities as time ed that the slope of the EDF heated by Ohmic disspation.
ith respect to the maximum ymmetrically.
from the low altitude subtudes.Fig. 6, panel b, shows t velocities are v d = 0 v the v the at t = 90 ν −1 0 .This have been accelerated from the tails begin to sprawl as a cts: enlarge the EDF symmetrielectrons having positive veight asymmetry between the presence of runaway electrons.The EDF are quite different along the box, but they are all non-Maxwellian.The runaway electrons can be seen more easily at higher altitudes.In Fig. 7, panel a, we present the total current density as a function of time for the same altitudes shown in Fig. 6.We can see that for every altitude, the current density increases very rapidly at the beginning but slows as time proceeds and eventually becomes constant.At higher altitudes, the first stage is longer as it takes more time for the electrons of the bottom altitude to reach the top of the box.The final current density is around 700 µA.m −2 for each altitude, which is similar to those reported by Noël et al. (2000).
Then, we are interested in the runaway current density.In order to differentiate between the current carried by thermal -Ohmic dissipation tends to enlarge the EDF symmetrically and, -The runaway effect acts on electrons having positive velocities, thus creating a slight asymmetry between the two tails of the EDF.The tail in the semi-space v>0 is larger than the one in v<0.This is apparent at t=90 ν −1 0 .
The distortion of the EDF due to the runaway effect can be clearly seen at the upper most altitude (Fig. 6 panel c).The drift velocity increases from v d =0 v the at t=0 ν −1 0 to v d =0.25 v the at t=90 ν −1 0 .Freely accelerated electrons from the lower altitudes create a high energy tail of the EDF.After 25 ν −1 0 , the distribution is clearly non-Maxwellian.For t>50 ν −1 0 , the tails are very asymmetric.This confirms the presence of runaway electrons.The EDF are quite different along the box, but they are all non-Maxwellian.The runaway electrons can be seen more easily at higher altitudes.
In Fig. 7, panel a, we present the total current density as a function of time for the same altitudes shown in Fig. 6.We can see that for every altitude, the current density increases very rapidly at the beginning but slows as time proceeds and eventually becomes constant.At higher altitudes, the first stage is longer as it takes more time for the electrons of the bottom altitude to reach the top of the box.The final current density is around 700 µA.m −2 for each altitude, which is similar to those reported by Noël et al. (2000).
Then, we are interested in the runaway current density.In order to differentiate between the current carried by thermal electrons and the current carried by runaway electrons we need to compute the runaway current density.This is done by using Fig. 1.In Fig. 1, to the left of the dotted line corresponding to v e <2.37 v the =0.7 eV the frictional force is dominant therefore the electric force can be neglected.For velocities larger than 2.37 v the corresponding to the region on the right of the dotted line, the electric force dominates.As a consequence, the electrons are accelerated for velocities larger than 2.37 v the .These electrons are called the runaway electrons and the current that they carry is the runaway current density.
We use this critical velocity (corresponding to 0.7 eV) to determine whether or not the electron is a runaway and in order to calculate a runaway current density.In Fig. 7. panel b, we can see that the runaway current density increases with time at each altitude: the electric field accelerates the electrons and so the runaway current density increases.We observe that the maximum runaway current density is about 100 µ A/m 2 .The ratio of runaway current density to the total current density is also of interest (see Fig. 7, panel c).In this plot, we represent the evolution of the ratio of the runaway current density to the total current density as a function of time for each altitude.We see that the ratio of runaway continues to increase particularly at high altitudes.The runaway electron can carry up to 20% of the total current density.

Summary and conclusion
The primary goal of this paper was to study electrons moving through a simplified ionospheric gas composed of ions and neutrals under the influence of a static electric field.To achieve this goal, a model representing the collisions between the charged particles and neutrals was developed.The model consists of two parts.The first part involves a kinetic description based on the Langevin equation, whose coefficients are determined using the Fokker-Planck equation.This part of the model, considering e/e and e/i collisions, gives the new electron velocities and the new EDF.
The second part is a Monte-Carlo method using a "null collision" approach.This part of the model deals with the elastic and inelastic electrons/neutrals collisions.These collisions cause the energy loss of the electrons.These collisions also tend to isotropize the EDF.
We have presented two examples.The first one without electric field shows that the relaxation of the EDF towards a Maxwellian is realized with respect to the theoretical relaxation time.It also provides a good illustration of the different impacts of electrons/electrons and electrons/ions collisions in relation to electrons/neutrals collisions.The second example is a simplified model of the auroral ionosphere between 200 km and 300 km altitude.A decreasing parallel electric field from the bottom of the box to the top is applied with a maximum value of 0.03 mV/m.This value has been chosen in order to reproduce the results published by Noël et al. (2000).It is shown that a significant distortion of the EDF due to the runaway effect occurs on a time scale of about 20 ν −1 0 , corresponding to 0.25 s.In other words, the EDF are non-Maxwellian all along the simulation box although the distortion of the EDF is more pronounced at higher altitudes.The maximum current density calculated in this run is 700 µA.m −2 which agrees with Noël et al. (2000).Runaway electrons can carry up to 20 % of the total current density.Our results suggest that the conclusions of the fluid models could be significantly altered by kinetic effects such as the runaway effect.
This paper focussed mainly on the collision operators used to model the ionosphere.We have used these operators in a simple 1-D model of the ionosphere.In the future, the electric field will be computed self-consistently.In addition, a 2-D model is under study in order to include a perpendicular electric field and the perpendicular motion of the charged particles.We also plan to simulate altitudes below 200 km where the current closure takes place.Then we will be able to include the magnetic field.Since the ion distribution is assumed to be stationary in our model, the time-evolution of the ion mean temperature, density and drift velocity need to be computed.We plan to do this using the fluid equations with the electron parameters (n e , T e ) calculated from our kinetic model.

Fig. 1 .
Fig. 1.The coefficient of friction due to e/e and e/i collisions are shown in blue and red, respectively.The coefficients are plotted as a function of normalized electron velocity.The total coefficient of friction is shown in black.The dashed-dotted line shows where e/e and e/i coefficients of friction are equal.The dotted line depicts the limit where the electric force is equal to the frictional force.

Fig. 2 .Fig. 3 .
Fig. 2. Characteristics of the electrons for different times when we apply only electron/electron and electron/ion collisions: (a) the electron distribution function in logarithmic scale versus normalized squared velocity.The different colors represent different times during the simulation (b) the total kinetic energy as a function of normalized time (c) the entropy as a function of normalized time (c).

Fig. 4 .
Fig. 4.Schematic illustration of the model for a plasma in the presence of an electric field aligned in the z-direction.The electric field decreases with height.The particles' velocities are 3-D even though the model is spatially 1-D.When a particle hits one of the boundaries, it is re-injected into the system according to the electron velocity distribution function f hi e (v x , v y , v z ).These boundary conditions are described in details in Sect.3.1.

Fig. 5 .
Fig. 5.A diagram of the computational scheme used in the paper.
onal scheme used in the paper.

Fig. 6 .
Fig. 6.Velocity distribution functions as a function of normalized squared velocity.The distribution functions corresponding to three different altitudes (from the bottom to the top panel respectively 2.5, 50 and 100 km) are shown for different times (0, 1, 10, 25, 50 and 90 ν −1 0 ) using different colors

Fig. 6 .
Fig. 6.Velocity distribution functions as a function of normalized squared velocity.The distribution functions corresponding to three different altitudes (from the bottom to the top panel respectively 2.5, 50 and 100 km) are shown for different times (0, 1, 10, 25, 50 and 90 ν −1 0 ) using different color.

Fig. 7 .
Fig.7.The left column of panels represents the total current density as a function of normalized time.The middle column represents the runaway current density as a function of normalized time.The right column of panels shows the ratio of the runaway current density over the total current density as a function of normalized time.For the three columns, the top panels correspond to 100 km altitude, the middle panels to 50 km and the bottom panels to 2.5 km.