A hybrid simulation model for a stable auroral arc

. We present a new type of hybrid simulation model, intended to simulate a single stable auroral arc in the latitude/altitude plane. The ionospheric ions are treated as particles, the electrons are assumed to follow a Boltzmann response and the magnetospheric ions are assumed to be so hot that they form a background population unaffected by the electric ﬁelds that arise. The system is driven by assumed parallel electron energisation causing a primary negative charge cloud and an associated potential structure to build up. The results show how a closed potential structure and density depletion of an auroral arc build up and how they decay after the driver is turned off. The model also produces upgoing energetic ion beams and predicts strong static perpendicular electric ﬁelds to be found in a relatively narrow altitude range ( ∼ 5000–11 000 km).


Introduction
Dynamic interactions and energy flow between the magnetosphere and the ionosphere are associated with auroral emissions, which can be either diffuse or discrete.Quiet discrete aurora consists of one or several auroral arcs.During disturbed times, discrete aurora can be very complicated but usually can still be described as being a superposition of several auroral arcs, which however now may undulate, wiggle, and wind up.Thus, a single quiet auroral arc is a basic building block of ionosphere-magnetosphere coupling whose physical origin we should try to understand before the more complicated, dynamical phenomena can be efficiently studied.It is known that auroral arcs are elongated, narrow forms that are associated with inverted-V type electron precipitation (Evans, 1974;Lin and Hoffman, 1979).The term inverted-V originally referred to the shape of the dependence Correspondence to: P. Janhunen (pekka.janhunen@fmi.fi) of the electron energy on latitude but, nowadays is usually used to describe any electron energy spectrum having a distinctive monoenergetic peak; the latitudinal dependence does not need to resemble that of an inverted V letter (Newell, 2000).
It is known with high certainty that there is a U-shaped potential structure (Carlqvist and Boström, 1970;Lundin and Eliasson, 1991;Carlson et al., 1998) above the optical arc in the acceleration region.Thus in this paper we will use the term "arc" to refer not only to the optical arc but also its associated potential structure and its associated plasma cavity, if any.Magnetospheric electrons are accelerated downward and ionospheric ions upward by the upward parallel electric field associated with the potential structure (Lundin and Eliasson, 1991;Carlson et al., 1998).Ionospheric electrons cannot rise above the bottom of the potential structure and magnetospheric ions are prohibited from reaching the ionosphere by the potential, unless they are energetic enough to overcome the potential.Typically a density depletion is observed above the bottom of the U-potential (Persoon et al., 1988;Carlson et al., 1998) and several types of plasma waves are commonly observed at different parts of the potential structure (Gurnett, 1991;André, 1997;Lysak and André, 2001).
The potential contours that are U-shaped at low altitude must close somewhere.In principle this can happen either in the other hemisphere in a symmetric potential structure or at some intermediate altitude below the equatorial plane (Hallinan and Stenbaek-Nielsen, 2001).If the potential contours close in the opposite hemisphere, electrons which are accelerated downward at the bottom of the potential must enter the flux tube by moving perpendicularly against an electrostatic force somewhere in the magnetosphere.On the other hand, if the contours are closed below the equatorial plane, there must be a downward electric field at that altitude which the electrons must overcome in order to enter the potential structure.Thus, if the question of potential closure can be resolved, the class of possible models for powering the auroral arcs by magnetospheric mechanisms can be restricted in an important way.Recent results support the idea that the potential contours close below the equatorial plane (Janhunen et al., 1999;Janhunen and Olsson, 2001;Hallinan and Stenbaek-Nielsen, 2001).
Most attempts to understand the formation of auroral arcs have been based on simulation models.Magnetohydrodynamic (MHD) based three-dimensional simulations coupled with models of precipitation-induced ionisation production in the ionosphere have been developed by several authors (Sato, 1978;Watanabe et al., 1993).Some of these works also contain an anomalous resistivity added (Otto and Birk, 1993;Watanabe et al., 1993).In the MHD simulations the main emphasis has been in studying the feedback instability, i.e. the nonlinear feedback between ionospheric conductance reacting to electron precipitation and the field-aligned current carried by the precipitating electrons.Magnetic reconnection that might occur in the acceleration region, because of the presence of anomalous resistivity leading into elongated current sheet formation, has been studied with nonideal MHD (Otto and Birk, 1993).Both the feedback instability and anomalous resistivity were found to produce elongated field-aligned current sheets that were interpreted as the large-scale manifestations of auroral arcs.MHD which contains both the ionospheric feedback instability and anomalous resistivity has been found to produce U-shaped potential structures under some conditions (Watanabe et al., 1993).A model, having two electron fluids and perpendicular ion motion with Alfvén waves, was used to show that a model based on Alfvén waves produces a time-asymptotic state that can be regarded as electrostatic (Rönnmark and Hamrin, 2000).Finally, two-dimensional MHD-type model augmented with electron inertia has been used to study the development and striation of auroral arcs (Seyler, 1990).
One-dimensional, fully kinetic, simulations for understanding the origin of the parallel electric fields have also been developed (Ergun et al., 2000;Schriver et al., 2001).It was found that the bottom of the U-potential may consist of two separate layers of parallel electric field, the lower one of which separates the cold ionospheric plasma from the auroral density depletion and the upper one separates magnetospheric protons from the auroral density depletion (Ergun et al., 2000).Depending on the magnetospheric source plasma boundary conditions, both upward and downward parallel electric fields were found to develop (Schriver et al., 2001).Recently, a spontaneous formation of several fieldaligned current sheets resembling auroral arcs was found to occur using a two-dimensional semiglobal hybrid simulation (Swift and Lin, 2001).
In this paper we present a new type of two-dimensional electrostatic hybrid simulation model in which both ionospheric and magnetospheric electrons are assumed to follow a Boltzmann response, ionospheric ions are treated explicitly as particles and magnetospheric ions are, for simplicity, assumed to be so hot that they are unaffected by the potential structure and therefore form a uniform background.This type of simulation is able to describe the formation of parallel electric fields, auroral density depletion and upflowing ion beams self-consistently in a two-dimensional setting.Waveparticle interactions of electrons are not self-consistently included in the model but effects describing them can be added in the manner described below.Our model contains the same physics as the one-dimensional models of Ergun et al. (2000) and Schriver et al. (2001), except for the omission of electron kinetic effects which, however, are not essential for the questions we study here.Because of our assumption that magnetospheric ions are very hot, we cannot expect to reproduce the double transition layer feature of Ergun et al. (2000).The advantages of our model are that it contains full ion kinetics and two spatial dimensions together with the possibility of studying the global effects of assumed electron wave-particle interactions.With our model we cannot, however, describe the magnetospheric auroral arc powering mechanism self-consistently as the MHD simulations do; energy input from the magnetosphere is contained in the assumed electron anisotropy caused by wave-induced parallel energisation.The development of a new type of ion pusher ("monopole solver") increases the efficiency of our model considerably and thus increases its practical applicability.We verify that the model produces the observational satellite signatures of an inverted-V region such as ion beams, convergent electric fields, parallel electric fields and an auroral density depletion.In this paper we use the model for studying the closure of the potential contours and the lifetime of auroral density depletions after the inverted-V precipitation that created them has disappeared.
The structure of the paper is such that we first present the simulation model, then show some results with figures with a discussion of their implications.At the end there is a discussion and summary section.

Simulation model
To simulate a system composed of auroral flux tubes we separate the electron and ion time scales.Electrons are assumed to move so fast compared to ions that the electron distribution reaches a quasi-stationary state before ions start moving.Thus the problem breaks down into (1) computing the time-asymptotic electron density profile when the ion density is known and fixed, (2) moving the ions one step further in known electric and magnetic fields.The magnetic field is assumed to be fixed and equal to the dipole field all the time.Thus the model is electrostatic and the electric field is a potential gradient, E = −∇V .
Because we are interested in scale sizes larger than the Debye length, quasi-neutrality holds, i.e. the total electron and ion densities are equal everywhere and at all times.Thus there are four particle populations: magnetospheric ions, magnetospheric electrons, ionospheric ions (taken to be O + ), and ionospheric electrons.We denote their densities by n e , respectively.The ionosphere is assumed to be a complete absorber for both electrons and ions and secondary particle production is neglected.We now discuss each of these populations in turn.
Magnetospheric ions (n (m) i ) are assumed to be so energetic that they are virtually unaffected by the electric field and thus they form a static background positive charge density profile n (m) i .In reality, this assumption is often valid and is made here for the sake of simplicity.We assume that the hot ions form an isotropic Maxwellian population at high altitude and that ions hitting the ionosphere are lost completely.The presence of ionospheric loss modifies the density profile from a constant value so that the density is given by where n (m) src is the magnetospheric source plasma density, B is the local magnetic field and B (i) is the ionospheric magnetic field (Janhunen, 1999, Eq. (31), henceforth called J99).
Similarly to magnetospheric ions, also magnetospheric electrons (n (m) e ) are assumed, in the ground state of the system, to originate from an isotropic Maxwellian source plasma at high altitude with complete loss occurring at the ionosphere.In the absence of an electric field the magnetospheric electron density profile n (m) e0 is, because of isotropy, equal to the magnetospheric ion density profile n (m) i (Eq.1).When an electric potential profile V is present, the magnetospheric electron density is modified by the Boltzmann factor so that n is the source plasma temperature of magnetospheric electrons.The formula is only valid if the potential structure is sufficiently simple so that particle accessibility effects do not arise.This will be the case for the potential structures that form in our simulation.
Ionospheric ions (n i ) are the ones that are treated as macroparticles in the simulation.Thus their density is obtained, at each time step, directly by the charge accumulation procedure.The ions are emitted from the bottom of the simulation box corresponding to an isotropic Maxwellian population whose density and temperature is given below.
Ionospheric electrons (n e ) are, similarly to the magnetospheric ones, assumed to follow the Boltzmann response, n The magnetospheric source plasma density in all runs is n (m) src = 0.3 cm −3 .There are 105 grid points along the magnetic field direction and 60 in the perpendicular.The grid spacing varies because of the dipole geometry but is about 350 km in the field-aligned direction and about 370 m in the latitudinal direction near the ionosphere.The timestep is 0.7 s.The total number of macroparticles varies but is typically about 580 000.The duration of one 1800 s run with these parameters is about 6 h (700 MHz AMD Athlon processor).The other parameters are listed in Table 1.
The quasi-neutrality condition now reads It is assumed here that the magnetospheric source plasma electrons are in the same potential as the ionospheric electrons.If this is not what is wanted, the equation should be correspondingly modified.For each field line, this equation must be solved numerically for V at each altitude.The right hand side of the equation is known since n (i) i is accumulated from the simulation macroparticles representing ionospheric ions (remembering to take into account also flux tube scaling).Since the left hand side of Eq. ( 2) is a monotonic function of V , the solution is unique.We have now completed solving the electron response problem, under the various assumptions stated; i.e. from a known n (i) i we know how to calculate V and thus the electric field E = −∇V everywhere.
The remaining task is to push the ions one step forward in known electric and magnetic fields.Gravitation is also taken into account and can be handled in the same way as the electric field.The electric field E is the potential field −∇V and the magnetic field B is the dipole field.This is, in principle, an easy problem for which many methods are available.One of the simplest methods is the Buneman solver (Hockney and Eastwood, 1988), which is the staggered leapfrog method applied to the Lorentz force expression.A property of the Buneman solver is that it gives the correct particle drift velocity even if the timestep is longer than the Larmor period; but the Larmor radius is then exaggerated because the particle moves along a straight line during each timestep.The exaggeration of the Larmor radius is harmful only if it exceeds the perpendicular grid spacing.Thus we obtain the timestep condition t < x ⊥ /v i for the Buneman solver.Here x ⊥ is the perpendicular grid spacing and v i is the (maximum) ion perpendicular speed.
At the ionospheric end of the simulation box the perpendicular grid spacing x ⊥ is so small (∼ 1 km) that using the Buneman solver would require using a short timestep, making the simulation slow.One can obtain a more efficient scheme by approximating the electric and magnetic fields in the vicinity of the particle in such a way that an exact solution can be obtained.It turns out that approximating B by the field of a hypothetical magnetic monopole gives a good method which is outlined in Appendix A. In the monopole solver the timestep condition is t < x /v i , where x is the parallel grid spacing.Thus one can use a factor of x / x ⊥ longer timestep than in the Buneman solver; this is a significant improvement because x / x ⊥ ∼ 10 3 near the ionosphere (see below).
When the ionospheric ions reach higher altitudes, in a stationary state the conservation of particles dictates nv/B=const where n is the density, v is the parallel veloc-ity and B is the magnetic field.This gives a useful rough estimate even though the state in the simulation is not completely stationary.Thus n ∼ B/v, so that both flux tube expansion and ion energisation contribute to a rapid lowering of the ionospheric ion density as a function of altitude.If nothing is done to improve the situation, most computing power is spent at low altitudes and the ion distribution function in the interesting high altitude regions is badly undersampled.Our solution to this problem is to split a particle into two whenever the grid cell, where the particle resides, contains too few macroparticles.Each macroparticle carries a number (weight) telling how many physical particles it represents.After each splitting, the weights of the daughter particles are half of the weight of the original unsplit particle (Kallio and Janhunen, 2001).The reverse procedure, joining particles together if a grid cell contains unnecessarily many macroparticles, is not needed here since the ions are almost always moving upward and thus along a widening flux tube.
We have now outlined the solution method of the coupled ion and electron equations but not yet discussed how the auroral arc is actually driven.To this end we take a rather general approach and assume the existence of a primary electron density n (m),extra e which is due to electron parallel energisation, presumably caused by wave-particle interactions.Thus we replace the quasi-neutrality condition, Eq. (2), by From observations it is known that the electron plasma on auroral field lines at ∼ 4 − 6 R E radial distance is often composed of a "middle-energy" population (100-1000 eV) and a hot population (a few keV) (Janhunen et al., 2001;Janhunen and Olsson, 2000).The middle-energy population is often anisotropic in such a way that the parallel temperature T is larger than the perpendicular T ⊥ while the hot population is isotropic.The middle-energy anisotropies are correlated with bursty plasma waves.In Eq. (3) we take the n (m) e0 term to represent the hot population and n (m),extra e the middle-energy population.Our working hypothesis is that a resonant electron parallel energisation mechanism affects the middle-energy population only, thus producing the T > T ⊥ anisotropy while leaving the hot population intact.
To get an explicit model for n (m),extra e we compute the electron density altitude profile for a bi-Maxwellian source plasma distribution function f BM , where W and W ⊥ are the parallel and perpendicular kinetic energies, respectively, T and T ⊥ are the temperatures, m e is the electron mass and n (m),extra e0 is the source plasma density associated with the anisotropic component at the magneto-spheric end of the simulation box.The distribution function f BM is normalised so that (5) (Eq.(3) of J99).To compute the density profile we utilise Eqs. ( 16) and ( 18) of J99 in the absence of a potential drop (V = V 0 = 0 in that paper's notation).Note that there is a factor 2π/m 2 missing in Eq. ( 16) of J99.The density is obtained by integrating over the parallel and perpendicular electron energies E and E ⊥ , Here b = B/B (m) and b 0 = B (i) /B (m) , where B (m) is the magnetic field at the magnetospheric end of the simulation box and B (i) is the ionospheric field.The upper limit of the parallel energy integration in Eq. ( 6) depends on b 0 ; parallel energies larger than the limit are in the loss cone and thus excluded from the integration.The parameter b 0 thus determines the loss cone width and by setting b 0 → ∞ we would obtain a case without ionospheric losses.The total density n (m),extra e is composed of upgoing and downgoing parts which are given by n up = n, n down = lim b 0 →∞ n (see J99).Doing the integrals we obtain, after some simplifications, One sees that whenever T = T ⊥ , Eq. ( 7) produces an altitude-dependent n (m),extra e , from which it follows that the potential V must also depend on altitude in a specific manner in order to satisfy Eq. (3).In the absence of anisotropy and without ionospheric losses (T = T ⊥ and b 0 → ∞), Eq. ( 7 grows downward because of flux tube convergence, except that in the vicinity of the ionosphere it gets reduced again because of the ionospheric loss.The reduction close to the ionosphere would be smaller if secondary electrons were taken into account. The strength of the arc formed can be controlled by the parameters T , T ⊥ and n (m),extra e0 . By increasing the T /T ⊥ ratio one makes the arc stronger, the potential structure deeper and its bottom altitude lower.A similar effect is also produced by increasing n (m),extra e0 . In the simulation here, we use Gaussian profiles for T /T ⊥ and n (m),extra e0 with half-width 0.05 o .This parameter defines the latitudinal width of the arc that develops.We assume that n (m),extra e is not affected back by the electric field.This is just a way of parameterising the driver, not an essential physical limitation.
The coordinate system used in the simulation is described in Appendix B. We use a coordinate system which is designed for allowing an exact dipole field without sacrificing numerical performance.
We initialise the runs by letting the simulation box first be completely empty and starting to emit ionospheric cold ions from the lower boundary.The emitted ions gradually move up and fill the simulation box with a density that drops approximately exponentially with altitude, since gravity lets only the highest energy ions reach high altitudes.During this initialisation phase we do not compute the electric field, i.e. the ions are treated as independent test particles.The length of the initialisation phase is 1 hour.We refer to t = 0 as the situation at the end of the initialisation phase, because at t = 0 the self-consistent run starts.
We now summarise conceptually the operations during one timestep in the simulation algorithm.The density due to ionospheric ions is first accumulated from the ion positions.This, summed with the density of hot magnetospheric ions (a fixed profile), gives the total plasma density which we require to be equal to the total electron density due to quasi-neutrality.Since the electrons follow a Boltzmann response, their densities depend on the potential V by Eq. (3).We solve V numerically from Eq. (3) to get the electric field E = −∇V , which is used to propagate the ions one step further.The magnetic field is constant and equal to the dipole field all the time.
At a point where the term n (m),extra e is nonzero, the electron density due to other electrons must be reduced in order to retain quasi-neutrality with the ion density at that point.The system accomplishes this reduction by creating a negative potential structure around the point which repels electrons and thus reduces their density locally.Qualitatively, one may think that introducing negative charge creates a neg- 3 Simulation results

Baseline run
The run described here is such that the driver is turned on at t = 100 s and turned off at t = 700 s, after which the run is continued for another 10 minutes to study the decay of the formed arc.
In Fig. 2 we show the density due to ionospheric ions at t = 400 s, i.e. 5 minutes after turning on the process producing the primary charge cloud.A density depletion has started to develop and has proceeded down to about 6500 km altitude.Ion outflow is also visible as an enhanced density in Fig. 2 at the centre of the arc up to about 30 000 km altitude.Since the ion beam ions move fast, their contribution to the density is small.The particle flux conservation requires that nv/B is constant where n is the ion beam density, v the ion beam velocity and B the local magnetic field.When the ion beam accelerates upward, v increases and B decreases so both effects combine to produce a rapid decrease of n.Notice that hot magnetospheric ions are not included in the density plotted here.The electric potential contours at t = 400 s are shown in Fig. 3.The lowest dip of the potential structure has proceeded down to 6500 km altitude at this time.The depth of the potential minimum is about 2.5 kV.The depth can be changed by changing the parameters of the process that causes the primary negative charge cloud (see previous section and Table 2).The potential structure appears somewhat more narrow than the density depletion.Most of the potential contours close below 20 000 km altitude.
The ionospheric ion density at t = 700 s (just before shutting off the driver) is shown in Fig. 4. At this stage, a clear density depletion starts from ∼ 5500 km upward.The density enhancement due to the ion beam is now clearly visible at the centre of the arc.The ions, that the potential structure has sucked up from the dense and cold ionospheric plasma layers, form the ion beam and have partly exited the system.Since the flux tube widens strongly as a function of altitude, the contribution of these ions to the plasma density at high altitude is small, however, and is invisible when the total plasma density (ionospheric+magnetospheric) is plotted (Fig. 5).Fig. 5 shows the lower portion of the density depletion, but above ∼ 15 000 km the total density is dominated by the magnetospheric ion background and thus does not show any structure.
The potential structure at 700 s, i.e. just before turning off the driver (Fig. 6) has moved downward by the same amount as the density depletion.The lowest dip reaches 5500 km altitude.The maximum potential depth is 3 kV.The perpendicular and parallel electric field components can be computed from the 10-s averaged potential and are shown in Figs.7 and 8, respectively.The perpendicular field reaches ∼ 400 mV/m magnitude in the acceleration region (∼ 7000 km altitude) and has a convergent signature.The parallel electric field is much smaller in magnitude (Fig. 8); the maximum upward parallel field of ∼ 1.6 mV/m occurs near the bottom of the acceleration region (∼ 6000 km altitude).Downward electric fields of smaller magnitudes occur in upper regions.Notice that strong perpendicular electric fields are concentrated in a rather narrow altitude range (∼ 5000-11 000 km).
After the driver has been turned off, the potential structure disappears in the electron time scale which is infinitely fast in the present simulation.However, the density depletion remains for ∼10 min.In Fig. 9 we show the density of ionospheric ions at t = 1300 s, i.e. 10 min after turning off the driver.Some remnants of the upgoing ion beam can still be seen around 20 000-30 000 km altitude.Below about 10 000 km, the plasma has almost returned to its original state but, above 10 000 km, a density depletion can still be identified.This stage of development could be called a "dead" auroral density depletion, as opposed to a "live" density depletion shown in the other figures.
In Fig. 10 we show simulated horizontal satellite passes through the system at t = 700 s (the fully developed potential structure, just before the driver is turned off).Satellite crossings are shown at 3000, 6000, 9000 and 12 000 km altitudes, including the total plasma density and the plasma potential along the orbit.At 3000 km altitude there is essentially no signal of the potential structure (the weak and wide potential minimum that is seen there is due to the ionospheric electric field).At 6000 km altitude one can see the formation of a sharp potential structure and a wider density depletion.At 9000 km the potential structure is wider, as is the density depletion.At 12 000 km the satellite is in the centre of the potential structure.

Other simulation runs
To investigate the effect of some input parameters in the baseline run we make four other runs.The runs and the changed parameters are listed in Table 1.The corresponding results for the potential bottom altitude, the depth of the potential, the density depletion minimum ionospheric ion density and the maximum perpendicular and parallel electric fields are given in Table 2.
In the "Anisotropy" run, the middle-energy electron anisotropy T /T ⊥ is lowered from 8 to 3. The result is a weak arc at very high altitude; the potential depth is only 1 kV and the bottom of the structure resides at 12 500 km altitude (Table 2).In a study of potential structures as a function of altitude using Polar data, only weak potential structures (less than 2 kV depth) were seen at high altitude (Janhunen et al., 1999).
In the "Lowdens" run the partial density of the middleenergy anisotropic electron population is halved.This brings about a rather similar outcome as the "Anisotropy" run, al- though now the potential structure bottom altitude is little lower (11 000 km) and the structure is more concentrated in altitude.
In the "Highdens" run the density of the anisotropic electrons is doubled.The resulting arc is a strong one, the potential depth being 12 kV and the bottom altitude as low as 2600 km.The peak perpendicular electric fields exceed 1 V/m.
In the "Winter" run, the ionospheric source plasma density is halved.The potential bottom altitude is now 2600 km but the potential depth (3.7 kV) is rather similar to the baseline run.The peak perpendicular field (800 mV/m) is twice as large as in the baseline run.
In all the runs presented here, no perpendicular potential difference was assumed to exist in the magnetospheric source plasma, i.e.Eq. ( 3) was used as is, without adding any offsets in the Boltzmann factors.We also made one run by assuming a Gaussian-shaped negative −3 kV potential profile in the source plasma, of the same half-width as the arc and without n (m),extra e .Thus, Eq. ( 3) was, in this run, replaced by  where V 0 is the Gaussian potential profile of the source plasma.The arc, in this case, is driven by magnetospheric shear flow instead of wave-induced parallel electron energisation and the resulting potential structure is U-shaped, i.e. the potential contours do not close within the simulation box.The low-altitude characteristics of the run are very similar to the baseline run but some differences occur at high altitude where the potential structures differ.In particular, convergent perpendicular electric fields appear not only at low altitude but also at high altitude in the shear-flow driven run.

Discussion and summary
Equation (3) determines the potential V when the densities of the various plasma components are known.To get some insight into the equation, let us consider it in two limiting cases and get an expression for the potential in each.Below the bottom of the acceleration region, the potential V is much smaller than the magnetospheric electron temperature, In this case the exponential multiplying n (m) e0 is almost unity and the equation can be solved for V to yield (remembering that n This expression gives a potential of the order of the cold ionospheric electron temperature T is no longer valid.When this occurs, however, we know that e so that the first exponential in Eq. ( 3) can be dropped and the equation can again be solved for V : . ( This produces potential values which are of the same order of magnitude as T (m) e and it can be used to estimate the depth of the potential structure.The parameter n (m),extra e appearing in these formulae depends on altitude (Eq.( 7) and Fig. 1) and so depends on the partial density of the middle-energy anisotropic component at the magnetospheric end of the simulation box n (m),extra e0 as well as on the middle-energy component temperature ratio T /T ⊥ .
We summarise our findings briefly.
1.The simulation is based on a separation of electron and ion time scales, as well as the separation of the parallel and perpendicular dynamics.
2. The "monopole solver" is needed to avoid resolving the ion travel time through a perpendicular grid cell (with the monopole solver, one only needs to resolve the travel time over a parallel cell).
3. After starting the driver, the potential structure and the associated density depletion form rather quickly but first reside at higher altitude.Within some minutes, the potential structure "digs" its way into the denser ionospheric plasma layers, making the ions flow up as fieldaligned energised beams.An almost stationary state is reached after about 10 minutes; at this time the lower boundary of the density depletion has reached so high density plasma that further downward motion is very slow.
4. In the baseline model, strong convergent perpendicular electric fields are confined in a relatively narrow altitude range (about 5000-11 000 km).
5. The strongest parallel fields (of about 1 mV/m magnitude) point upward.There are much weaker downward fields at higher altitudes, however.
6.After turning off the driver, the potential structure disappears immediately (in the electron time scale) but the density depletion and partly the ion beam remain ∼10 min.
7. The arc can be made stronger by increasing the anisotropy T /T ⊥ or the density of the anisotropic component n (m),extra e0 . A strong arc has a deeper potential structure and it extends to lower altitude than a weak arc.
8. The arc bottom altitude can be lowered by decreasing the ionospheric source plasma density n (i) src (winter-like conditions).Changing this parameter does not appreciably modify the potential structure depth.9.By adding an offset in the Boltzmann factor, the arc can also be driven by magnetospheric shear flow.In this case the potential contours do not close in the simulation box.The low-altitude phenomenology of the shear flow driven arc is very similar to the anisotropy driven arc of the baseline run.
The model described here is in accordance with the "cooperative" model of auroral acceleration presented earlier (Janhunen and Olsson, 2000).In this model, wave-particle interactions and a potential structure are together responsible for producing inverted-V electron signatures; waves do the energisation and the potential structure gives the shape of the low-altitude energy spectrum.In the present simulation, the wave-electron interactions are modelled only by their assumed effect on the middle-energy electrons (the anisotropic population n (m),extra e in Eq. 3).In the future, more realistic wave-particle interaction models should be incorporated in the simulation.A straightforward way of doing this is to model electrons as particles and solve the electric field from Poisson's equation.This requires a lot of computing power which, however, may become available in the rather near future.Another way to extend the current model would be to make it three-dimensional or to add the return current region adjacent to the arc by modelling the relevant ionospheric physics.

Appendix A: Monopole solver
The ion pusher needs to update the ion position and velocities in a known electric and magnetic fields.The electric field E is electrostatic, E = −∇ , and the magnetic field B is the dipole field.The convergent nature of the magnetic field should be taken into account as accurately as possible.We shall approximate the dipole field in the vicinity of the ion with the field of a magnetic monopole.The monopole field is convergent and is thus a much better approximation for the dipole field than a constant field.
Consider a charged particle in a magnetic monopole field and a radial electric field (Clemmow and Dougherty, 1969), where the hat denotes a unit vector, v = dr/dt, and the dot denotes a time derivative, ṙ ≡ dr/dt = v.
Forming the vector product of Eq. ( 11) with r and taking the time derivative, we obtain On the other hand one can derive Substituting Eq. ( 13) into Eq.( 12) we obtain We can integrate Eq. ( 14) once with respect to time to find that where k is a constant vector.
Taking the dot product of Eq. ( 15) with r we obtain i.e. r • k = 1/k.Since r • k = cos θ where θ is the angle between k and the position vector r, we infer that θ is a constant of motion and the motion of the particle therefore is confined to a conical surface (Clemmow and Dougherty, 1969).The symmetry axis of the cone is given by k.
Introducing a Cartesian coordinate system (x, y, z) with ẑ = k and the associated polar coordinates (r, θ, ϕ), θ does not change while the particle moves, only r and ϕ do so.Thus, to solve the orbit, it is enough to solve r(t),ϕ(t).Taking the cross product of Eq. ( 15) with r we find from which we obtain On the other hand, since the radius of the cone is r sin θ and because |r × k| = | sin θ | we can write The minus sign follows from the fact that, in a right-handed spherical coordinate system, ϕ grows in the direction of z×r.Combining Eq. ( 18) with Eq. ( 19) we obtain Using Eq. ( 20) we can solve ϕ(t) by integration, if r(t) is known.
To find r(t), let us form the quantity d(v 2 )/dt and utilise Eq. ( 11) and Eq. ( 13): from which we obtain, by integrating once, where U is a constant (the "total energy").Let us also form (d/dt)(r ṙ): where we used Eq. ( 22) in the last phase.Equation ( 23) is an ordinary second order differential equation for r(t) which is, however, difficult to solve analytically unless (r) has some special properties.For our purposes here we only need to find an analytical form of (r) which allows for explicit solution and contains enough free parameters to allow a local approximation which is valid in the vicinity of the particle.
To this end we consider the case when (r) contains terms of the form r ±2 only, i.e.
(r) = 1 2 where p 1 and p 2 are real parameters.Using this form for (r) we derive, from Eq. ( 23), the equation Notice that p 1 drops out from Eq. ( 25).Introducing a new variable u ≡ r 2 , Eq. ( 25) takes the form whose general solution is (returning again to using r(t)) where α and β are constants of integration.These constants must be determined from the initial conditions which gives us Substituting Eq. ( 27) into Eq.( 20) and integrating with respect to t, one can derive the orbit (r(t), ϕ(t)).If one wants to avoid complex numbers, one must consider the cases p 2 > 0 and p 2 < 0 separately.The integration is straightforward and yields arc-tangent functions in both cases.In the ion pusher one also has to differentiate the orbit in order to obtain the velocities v r (t) and v ϕ (t).The resulting expressions are not terribly long but, since they are not particularly illuminating either we shall not reproduce them in this Appendix.
One also has to find the position of the monopole as well as the parameters q, p 1 and p 2 that give a best fit to the wanted magnetic and parallel electric fields at a point.It turns out that one can reproduce not only B and E r but also (B • ∇)B and ∂E r /∂r in the vicinity of the point.We give the resulting expressions without derivation (the derivation is simple, though): where R is the distance from the point to the monopole (the direction of the monopole is equal to ± B), In these formulae E R is the parallel electric field.In order to account for the perpendicular electric field, we make a coordinate transformation to a frame moving with the E × B velocity with respect to the inertial frame, calculate the monopole orbit there, and transform the result back to the inertial frame.
We have verified that the monopole solver described in this Appendix produces accurate results by comparing the results with the Buneman solver (Hockney and Eastwood, 1988) with a very small timestep.We have also verified the monopole solver result with the exact solution in some special cases.6 Appendix B: Coordinate system Our coordinate system (χ, θ i , ϕ i ) is defined so that θ i and ϕ i are, respectively, the magnetic dipole colatitude and longitude of the ionospheric footpoint and χ is a coordinate along the magnetic field such that χ = 0 at the ionospheric plane.
The exact nature of χ will be defined below.The benefit of the (χ, θ i , ϕ i ) coordinate system is that the coordinates are "field-aligned" in the sense that a stack of grid cells having the same θ i and ϕ i belong to the same dipole flux tube.
Let us define a transformation from Cartesian coordinates (x, y, z) to (χ , θ i , ϕ i ) by going through ordinary spherical coordinates (r, θ, ϕ), defined by r 2 = x 2 + y 2 + z 2 , cos θ = z/r, tan ϕ = y/x.The equation of a dipole field line is r = L sin 2 θ , where L is the McIlwain L-parameter (the radial distance at which the field line intersects the equatorial plane).The L-parameter is related to the footpoint colatitude θ i by 1 = L sin 2 θ i .Thus sin θ i = sin θ/ √ r.We also have ϕ i = ϕ.
Let us denote the field-aligned distance from the ionosphere by s, s = ds.Defining χ = s would be a natural choice and would give χ a simple physical interpretation.However, it would make the coordinate system quite oblique at a larger distance from the Earth (for example, smaller L-values would reach the equatorial plane at much smaller s-values than higher L-values).Our next trial is to define χ = s/L.We can now calculate s as The function F (x) has an explicit expression, but we shall not use it in the simulation because it is rather slow to compute.The function F is plotted in Fig. 11.Thus, the mapping to (χ , θ i , ϕ i ) coordinates is ϕ i = atan(y/x), θ i = asin(sin θ/ √ r), and χ = F (cos θ i ) − F (cos θ) if χ is defined as χ = s/L.
In the ϕ i = 0 sector, the resulting (χ, θ i ) coordinate grid is shown in Fig. 12 when the magnetic latitude of the footpoint 90 o − θ i varies from 65 o to 73 o .We see that the coordinate grid is field-aligned as it should, but still becomes quite oblique at larger distances from the Earth.Although not a problem in principle, it is slightly inconvenient.
ionospheric electron density profile in the V = 0 case and T (i) e is the temperature of the ionospheric electrons.There are several possibilities for specifying a model for n (i) e0 , all of which, however, yield almost identical final simulation results.This is because constant value which is equal to the ionospheric source density n (i) src at the bottom of the simulation box.
e. a constant.Fig. 1 shows n (m),extra e for two different values of T /T ⊥ .Basically, n (m),extra e

Fig. 1 .
Fig. 1.The quantity n e (m),extra given by Eq. (7) as a function of radial distance R (in units of Earth radius) for T /T ⊥ = 2 (solid line) and T /T ⊥ = 8 (dashed line).Normalisation is such that n e (m),extra = 1 at R = 5.

Fig. 2 .
Fig. 2. Density of cold ionospheric ions at t = 400 s, averaged over the previous 10 s.

Fig. 9 .
Fig. 9. Density of cold ionospheric ions at t = 1300 s, averaged over the previous 10 s.

Fig. 10 .
Fig.10.Simulated satellite passes at 3000, 6000, 9000 and 12 000 km altitudes, showing the total plasma density and the plasma potential along the satellite orbit as a function of invariant latitude.
i defines the bottom of the acceleration region, above which the approximation |eV | T (m) e

Table 1 .
Run input parameters

Table 2 .
Run state at t = 700 s Run Pot.depth Pot.bottom alt.Minimum ionosph.dens.Max.E ⊥ Max.E