Vlasov Simulations of Trapping and Loss of Auroral Electrons

The plasma on an auroral field line is simulated using a Vlasov model. In the initial state, the acceleration region extends from one to three Earth radii in altitude with about half of the acceleration voltage concentrated in a stationary double layer at the bottom of this region. A population of electrons is trapped between the double layer and their magnetic mirror points at lower altitudes. A simulation study is carried out to examine the effects of fluctuations in the total accelerating voltage, which may be due to changes in the generator or the load of the auroral current circuit. The electron distribution function on the high potential side of the double layer changes significantly depending on whether the perturbation is toward higher or lower voltages, and therefore measurements of electron distribution functions provide information about the recent history of the voltage. Electron phase space holes are seen as a result of the induced fluctuations. Most of the voltage perturbation is assumed by the double layer. Hysteresis effects in the position of the double layer are observed when the voltage first is lowered and then brought back to its initial value.


Introduction
Auroral optical emissions are caused by electrons that have been accelerated by electric fields that are parallel to the magnetic field.These parallel fields have been measured, for example, by the S3-3 (Mozer et al., 1977), Viking (Lindqvist and Marklund, 1990), Polar (Mozer and Kletzing, 1998), and FAST (Ergun et al., 2002) satellites.The parallel elec-tric fields can be concentrated in thin electric double layers as was proposed by Alfvén (1958).Their presence has been confirmed by satellite-based measurements using the S3-3 (Temerin et al., 1982), Polar (Mozer and Kletzing, 1998), and FAST (Ergun et al., 2002) spacecraft.The latter also found strong double layers in the downward current region (Andersson et al., 2002).FAST data has also been used to study electron phase space holes and electromagnetic radiation in the vicinity of double layers (Pottelette et al., 2014).Extensive studies of parallel electric fields and their role in auroral acceleration have been performed using the Cluster spacecraft (Marklund et al., 2001(Marklund et al., , 2011;;Sadeghi et al., 2011;Li et al., 2013).Forsyth et al. (2012) used Cluster data to study a case with a time varying acceleration voltage and found that the majority of the change, V , was assumed below the spacecraft altitude of about 4700 km.Double layers have also been observed in laboratory experiments with magnetised plasmas (Torvén and Andersson, 1979;Torvén, 1982;Schrittwieser et al., 1992).While in a double layer the electric force is balanced by electron inertia, in an inhomogeneous magnetic field also the magnetic mirror force can balance the force from the parallel electric field (Alfvén and Fälthammar, 1963).Laboratory experiments on parallel electric fields in magnetic mirror fields have been carried out in Q-machines (Sato et al., 1986(Sato et al., , 1988) ) and modelled using particle in cell simulations (Ishiguro et al., 1995).It was found in experiments, simulations, and theory of discharge plasmas that electrons accelerated in an electric double layer can cause oscillating electric field spikes when entering a density gradient on its high potential side (Gunell et al., 1996a, b;Gunell and Löfgren, 1997;Löfgren and Gunell, 1998;Brenning et al., 2006).
Published by Copernicus Publications on behalf of the European Geosciences Union.Song et al. (1992) derived a stability criterion for double layers in a converging magnetic field, showing that a stable double layer position is possible when electrons are accelerated into a stronger magnetic field.Ergun et al. (2000b) used a static model to find solutions to the Vlasov-Poisson system over a distance of a few Earth radii along an auroral flux tube.The Vlasov equation was used in plasma physics by Vlasov (1938), although an equation of the same form had been used earlier to study the motion of stars (Jeans, 1915).Poisson's equation was published by Poisson (1813).Hwang et al. (2009) studied ion heating and outflow using a Vlasov simulation model of the downward current region.Other ways of modelling include drift-kinetic simulations that have shown electron heating and parallel electric field related to shear Alfvén waves (Watt et al., 2004).Vlasov simulations of double layer formation in unmagnetised plasmas were performed by Singh (1980), who also found electron phase space holes on the high potential side of the double layer (Singh, 2000).

H. Gunell et al.: Auroral electrons trapped and lost
It was pointed out by Persson (1966) that a unique potential profile V (z) cannot always be found.Whipple (1977) introduced an effective potential U (z) = qV (z) + µB(z), where q is the particle charge, µ its magnetic moment and B the magnetic flux density.If there is a maximum in U (z) between the source of the particle population and the point in space under consideration, call it z = z 1 , the distribution function cannot be determined uniquely by the local V (z 1 ) and B(z 1 ).With a maximum in U (z), knowledge of V (z) and B(z) for all points between z 1 and the source is required.Chiu and Schulz (1978) showed that the condition for there being no maximum is that dV /dB > 0 and d 2 V /dB 2 ≤ 0, and they went on to consider only steady state solutions that fulfil that condition.If, instead of a maximum, a minimum exists in U (z), this constitutes an effective potential well, in which particles can be trapped.The trapped particle populations inhabit a region of phase space from which they cannot reach the system boundaries for the given shape of U (z). Equivalently, this region of phase space cannot be reached by particles entering through the system boundaries.Thus, if a steady state solution is sought, assumptions must be made about the trapped particle populations.Such assumptions were made in previous studies by, for example, Chiu and Schulz (1978) and Ergun et al. (2000b).The influence of the trapped electrons on V (z) was studied by the latter and by Boström (2004) by varying the density assumed for the trapped population.
Since auroral potential differences are changing with time, we allow the potential to evolve, and particles can be caught and trapped during the formation of a potential profile.Gunell et al. (2013a) found that an electron population became trapped between a double layer and the magnetic mirror field at lower altitudes.It has also been found that fluctuations are an important influence on the electron distributions, leading to electron conics such as those observed by the Viking satellite (André and Eliasson, 1992;Eliasson et al., 1996).Temporal variations of the potential have been observed by the Cluster spacecraft (Marklund et al., 2001(Marklund et al., , 2011;;Sadeghi et al., 2011;Forsyth et al., 2012).In this paper, we perform numerical experiments on a system that contains trapped electrons in order to study how they are trapped and how they can be released again.

Simulation model
We study the trapped and precipitating electron populations through the use of a Vlasov simulation code that is onedimensional in space and two-dimensional in velocity space (Gunell et al., 2013a).This model was used to model the plasma on an auroral magnetic field line and also to investigate how auroral acceleration can be simulated in a laboratory experiment (Gunell et al., 2013b).
In this model, the distribution function is f (z, v z , µ, t), where z is the spatial coordinate parallel to the magnetic field, v z is the parallel velocity, and µ = mv 2 ⊥ /2B(z) is the magnetic moment.The magnetic moment is a constant of motion, and therefore we have μ = 0.The Vlasov equation (Vlasov, 1938) that needs to be solved for our system is (Gunell et al., 2013a) In the gravitational acceleration a g only the component parallel to the magnetic field is included.The electric field is found by solving a Poisson type equation (Poisson, 1813) that has been adapted to the geometry of the problem at hand.The equation used here is (Gunell et al., 2013a) d dz where S is the cross section of the flux tube at the point where B = B S .The line charge ρ l on the right hand side of Eq. ( 2) is the charge per unit length of the flux tube.It is computed by integrating the distribution function over velocity space and forming the sum over all species s (Gunell et al., 2013a) In Eq. ( 2), the constant = 0 r includes an artificial relative dielectric constant r .This has been introduced for computational reasons; because λ D ∼ √ r and ω p ∼ 1/ √ r it allows the simulation to be run with a longer time step on a coarser spatial grid.Thus the system can be filled with plasma quickly at a high r value, and a time accurate simulation can subsequently be run with a lower r .The numerical experiments presented in this article were run at r = 4.97 × 10 4 .
For relativistic particles, instead of conservation of µ, the conserved quantity is (Alfvén and Fälthammar, 1963) where m 0 is the particle mass at rest, and γ is the Lorentz factor.In the simulations described in the following sections, the electrons that have their source at the magnetospheric end of the system are treated relativistically, and we therefore write the distribution as a function of µ rather than µ for this species.For further details about the model we refer the reader to our previous paper (Gunell et al., 2013a).

Initial state
Using a dipole model of the magnetic field, we have simulated the plasma on a flux tube of the L = 7 shell from the equatorial magnetosphere, where z = 0, down to the ionosphere, where z = 5.5 × 10 7 m.At the magnetospheric boundary of the system we impose a Maxwellian plasma with k B T e /e = 500 V and k B T i /e = 2500 V.At the ionospheric boundary we assume Maxwellian distributions with k B T /e = 1 V for both electrons and ions.All the ions included here are protons.We treat particles that have their sources at the two ends of the system as different species, and thus we have four species in the simulation: electrons from the magnetosphere; protons from the magnetosphere; electrons from the ionosphere; and protons from the ionosphere.The initial state of the numerical experiments that are reported in Sect. 4 is a simulation run with a voltage of 3 kV over the system.The initial state is illustrated in Fig. 1; Fig. 1a shows the plasma potential, and it is seen that about half of the voltage falls in a double layer at z = 5 × 10 7 m.The position can be determined as the position of the large negative spike in Fig. 1b.While half of the total voltage is concentrated in a double layer, the other half is in an outstretched region above it.Similar weak parallel electric fields have been found in simulations of Alfvén waves (Vedin and Rönnmark, 2008;Watt and Rankin, 2009).Our simulation is completely electrostatic and can therefore not simulate Alfvén waves.However, the electron time scale is faster than the time scale of Alfvén waves.Thus, an electrostatic equilibrium may be established in a few seconds, and that equilibrium be perturbed by Alfvén waves on time scales of tens of seconds or minutes.
The simulation leading to the initial state in this paper differs from the final state of the simulation by Gunell et al. (2013a) in that the µ resolution of the magnetospheric electrons has been improved, as can be seen in Fig. 1g, which shows the distribution function f (µ , v z ) near the double layer on its high potential side at z = 5.1 × 10 7 m.Here, and in all the following figures that show f (µ , v z ), the white lines indicate the boundaries of different phase space re- (g) Phase space density f e,M (µ , v z ) for the electrons from the magnetosphere at z = 5.10 × 10 7 m, which is close to the double layer on its high potential side.The phase space regions separated by the white lines are I: precipitating electrons; Ib: up-going electrons that will reflect and then precipitate; II: electrons which can reach the equator; and III: trapped electrons.The unit for f e,M (µ , v z ) is m −6 A −1 s.In panels (c)-(f), the colour scales have been normalised so that integrals over all v z give n s /B.The unit for f (z, v z ) is m −4 T −1 s.
after first being reflected by the electric field.In region II, the electrons can reach the equator, and region III contains electrons that are trapped between the magnetic mirror and the potential drop.The lines are calculated according to Gunell et al. (2013a), and take into account V and B at all z points between the point shown, z = 5.1 × 10 7 m in this case, and the boundary in question.Given enough time, electrons can reach the ionosphere if and they can reach the magnetospheric end of the system if where z = z 1 is the point in space under consideration, V (z 1 ) = V 1 , B(z 1 ) = B 1 , and z = L z is the z coordinate of the ionosphere (Gunell et al., 2013a).A reduced distribution function f (z, v z ), where we have integrated over the µ dimension, is shown in Fig. 1c-f for each species.
The relevant time scales can be estimated by computing approximate transit times for a few typical particles.For a magnetospheric electron with an energy of 500 eV it takes 3.8 s to cover the distance of 5 × 10 7 m to the double layer, and having been accelerated to 3 keV in the double layer, a precipitating electron traverses the high potential side to the ionosphere in 0.15 s.It takes a 2.5 keV ion 72 s to move from the magnetospheric end of the system to the double layer.Ions of ionospheric origin that are accelerated to 3 keV in the double layer reach z = 0 in 66 s.

Numerical experiments
Starting from the initial state described in Sect.3, we have performed three numerical experiments where the total voltage imposed on the system was changed during the simulation run.The voltage as a function of time for the three experiments is illustrated in Fig. 2.
In experiment 1 the voltage was first decreased from 3 kV to 1 kV during one second, kept constant at 1 kV for two seconds, increased back to 3 kV during one second, and then the simulation was run at 3 kV for another three seconds until t = 7 s.Experiment 2 is the opposite of experiment 1, that is to say, the voltage was first increased to 5 kV and subsequently brought back down to 3 kV.In experiment 3 the voltage was brought down all the way to 0 V before it was increased again to 3 kV.

Experiment 1
Figure 3 shows the phase space density f e,M (µ , v z ) in experiment number 1 at z = 5.1 × 10 7 m, which is a position close to the double layer on its high potential -low altitude (2) (3) -side.Panels (a)-(h) show one distribution function per second, starting with the distribution in the initial state.To the right of each panel is a time profile of the total voltage over the field line with the voltage at the time for which the distribution is shown marked "×".At t = 0 there was one population of precipitating electrons and one of trapped electrons.These trapped electrons were only present for µ 8 × 10 −12 A m 2 .The voltage was ramped down linearly from 3 kV to 1 kV, which was the voltage at t = 1 s.It is seen in Fig. 3b that the phase space region populated by electrons had shrunk somewhat at t = 1 s compared to the state at t = 0, and that the precipitating electrons had become slightly slower in v z .However, the more obvious change is that the boundaries between the different phase space regions, shown by the white lines, moved.At the lower voltage the region of phase space where electrons could reach the ionosphere was smaller, the region where electrons could reach the magnetosphere larger, and the part of phase space that contained the trapped electrons had undergone the greatest diminution of them all.
The electrons that at t = 1 s had found themselves in the region where they were allowed to reach the magnetosphere were lost during the next two seconds while the voltage was kept at 1 kV, and it can be seen that the phase space density in this region decreased for t = 2 s and t = 3 s, in panels (c) and (d) respectively, compared to the initial state.It is also seen that electrons had entered the previously empty region that can contain trapped electrons at low µ .
Between t = 3 s and t = 4 s the voltage was increased back to 3 kV, and then it was held at that 3 kV for the remainder of the simulation run.The boundaries between the different phase space regions returned approximately to their initial positions, but it can also be seen that the white curve reaches farther toward higher µ values in panel (h) than in panel (a).The region of phase space from which electrons had been lost when the voltage was lowered was refilled when the voltage was raised again between t = 3 s and t = 4 s.The highest phase space density of the initial state in the approximately triangular region between µ ≈ 8 × 10 −12 A m 2 and µ ≈ 4 × 10 −11 A m 2 was, however, not recovered completely.There was some increase in the phase space density for µ 2 × 10 −11 A m 2 when the voltage was constant in the interval 4 s ≤ t ≤ 7 s due to electrons moving in the z dimension.As µ is a constant of motion, there can be no flux in the µ dimension; a flux in the v z dimension, corresponding to electron acceleration, would be possible in principle, but no signs of such an acceleration can be seen in Fig. 3e-h.The electrons that appeared at low µ values when the voltage was lowered remained there throughout the run of the simulation.
To illustrate the electron motion along the magnetic field, Fig. 4 shows distribution functions f e,M (z, v z ) for electrons of magnetospheric origin.The left-hand column (Fig. 4a) shows a distribution that has been integrated over µ < 8 × 10 −12 A m 2 , and the mid column (Fig. 4b) shows the integrated distribution for µ ≥ 8×10 −12 A m 2 .The colour scales have been normalised so that integrals over all v z yield partial density divided by magnetic field, that is to say, where n p (z) denotes the partial density of magnetospheric electrons for the µ range in question.The vast majority of the electrons represented by the initial distribution shown in Fig. 4a 0 were precipitating.Only a small fraction was reflected by the magnetic mirror field and accelerated toward the double layer, where the electrons were slowed down.Then they were further decelerated by the weaker electric field at higher altitude until they returned downward to the double layer again.Thus, they were trapped between the magnetic mirror and the potential structure, specifically the part of the potential structure which was situated in an extended weak field region above the double layer.For the higher µ range shown in Fig. 4b 0 there were also electrons at lower |v z | that were trapped between the magnetic mirror and the double layer electric field.This can be seen by the sharp edge at the double layer position.
As the voltage was lowered, the potential drop located in the double layer decreased, but it did not disappear completely.This can be seen in Fig. 5a, where the green curve shows the initial potential profile and the red curve the potential at t = 1 s in experiment 1.As a result we see that part of the trapped population had been lost at t = 2 s, shown in Fig. 4b 1 , and as it takes some time for the electrons to turn around and leave the high potential side of the double layer, this population continued during the interval 1 s ≤ t ≤ 3 s, when the voltage was constant at 1 kV in Fig. 4b 2−3 .In the low µ range, electrons were reflected by the magnetic mirror, and as the voltage had been decreased they were not decelerated as much in the double layer as in the initial state.Therefore there is a gap between the upgoing and downgoing populations in Fig. 4a 2 for z 5×10 7 m.There was also a part of the low µ electron population that did not pass back through the double layer toward lower z values and instead got trapped on the high potential side.These were the same particles as those that are seen filling the trapped particle region for µ 8.0 × 10 −12 A m 2 in Fig. 3c-h.
In a steady state, trapped particles cannot reach the system boundaries, and, conversely, particles at the boundaries cannot reach the regions of phase space containing trapped particles.When the system is not in a steady state, particles can both enter and leave the traps.The empty part of the trapped region in Fig. 3b below µ ≈ 8 × 10 −12 A m 2 was no longer empty one second later in Fig. 3c, in spite of the region boundaries being largely the same in the two figures.Fig. 6a shows the phase space density f e,M (µ , v z ) at z = 5.1×10 7 m for t = 1.5 s -that is to say, between the times shown in Fig. 3b and c -and here the accelerated electrons that were in the precipitating region for the distributions shown in Fig. 3 had been decelerated into the trap.The reason for this is the fluctuations that appeared as a result of the imposed modulation of the voltage.Figure 6c shows the total voltage as a function of time (blue curve) together with the plasma potentials at two points near the double layer that initially both were on its low potential side, namely z = 4.92 × 10 7 m (black curve), z = 4.975 × 10 7 m (red curve).When the voltage first was lowered, all three curves decreased, but there was a rebound for z = 4.92 × 10 7 m so that at t = 1 s, the potential there was higher than at z = 4.975 × 10 7 m.This can also be seen by the dip in the red curve in Fig. 5a between the two vertical dashed lines in that figure.The oscillation that ensued peaked at a potential near 2 kV for z = 4.92 × 10 7 m shortly before t = 1.5 s.After that, the oscillation decayed, and the contribution from the increase of the voltage between t = 3 and 4 s to this oscillation was smaller than that of the initial voltage decrease.The red curve showing the potential at z = 4.975 × 10 7 m attaches itself to the blue curve around t = 2 s, showing that the double layer moved toward lower z values so that z = 4.975 × 10 7 m was on the high potential side after t = 2 s.The peak of the black curve corresponds to an electric field that was directed in the positive z direction, and that is what decelerated the electrons in Fig. 6a.
Figure 6b shows the partial phase space density f e,M (z, v z ) of magnetospheric electrons, integrated over 0 ≤ µ < 8.0 × 10 −12 A m 2 for t = 1.5 s.On the right hand side electrons being reflected by the magnetic mirror field are seen, and in the centre of the image electrons were reflected by the reversed electric field that existed in this region during the local potential peak shown by the black curve in Fig. 6c.This created two beam-like populations moving in the negative z direction.These beams can be seen at t = 2 s as the somewhat complicated structure at negative v z values in Fig. 4a 2 .

Experiment 2
In experiment 2 the voltage was first increased from 3 kV to 5 kV, held constant there for 2 s and then decreased back to 3 kV, where it was maintained until t = 7 s. Figure 7 shows the phase space density f e,M (µ , v z ) in experiment 2 at z = 5.1 × 10 7 m, which is a position close to the double layer on its high potential side.This is the same position that was shown for experiment 1 in Fig. 3.
As the voltage was increased, the boundary of the trapped electron region expanded both toward larger µ and larger |v z | values.The phase space density is seen to have increased just inside the parabola in Fig. 7b, both at its tip on the right-hand side of the image and for 1×10 −11 A m 2 µ 3× 10 −11 A m 2 .The electrons that were trapped by the expansion of the trapped particle region had entered the high potential side of the double layer, and after being reflected by the magnetic mirror they encountered a higher electrostatic potential barrier due to the increased voltage -as is shown by the plasma potential at t = 1 s in Fig. 5a (black curve) -which prevented their return to the low potential side.When the voltage subsequently was returned to 3 kV again the newly trapped electrons were released from the trap.
In  The initial distribution for µ < 8.0 × 10 −12 A m 2 in Fig. 8a 0 is the same as that in Fig. 4a 0 , as both experiments started from the same initial state.Most of the electrons were precipitating, and the small fraction that was not can be seen to be trapped between the magnetic mirror and the potential structure above the double layer.After the voltage was increased, the part of this small trapped distribution that was located above the double layer was lost, as an increased acceleration in the double layer allowed them to precipitate.The black curve in Fig. 5a shows the plasma potential at t = 1 s in experiment 2, and it is seen that the double layer voltage increased when the total voltage was raised.After the voltage was returned to 3 kV the region that originally contained the trapped electrons was filled anew.The distribution between z = 4 × 10 7 m and z = 5 × 10 7 m was more structured at t = 7 s than the original distribution, and that is a result of the potential structure that appeared in that region after the return of the total voltage to 3 kV.Similar structuring is also seen for the intermediate and high µ ranges.A fraction of the low µ electron distribution can also be seen returning toward the magnetospheric side of the system in Fig. 8a 6 and a 7 .This was caused by the decreasing double layer voltage during the lowering of the total voltage between t = 3 s and t = 4 s.Electrons were accelerated towards higher z by a high double layer voltage, then they were reflected by the magnetic mirror, and when they reached the double layer again its voltage had decreased, enabling them to reach the magnetosphere.
This last phenomenon can also be seen in Fig. 8b 5−7 for the intermediate µ range.The initial distribution in this range was symmetrical with respect to v z = 0 (Fig. 8b 0 ).At the higher voltage in Fig. 8b 1−3 the negative v z part was diminished as the higher voltage allowed more electrons to precipitate rather than be reflected.The initial symmetry did not return completely during the simulation run.When the voltage had been returned to 3 kV at t = 4 s the changed potential profile caused structuring and there were fluctuations causing heating of the electrons moving in the negative z direction.The double-peaked distribution for intermediate µ that was noticed in Fig. 7f-h also appears in Fig. 8b 5−7 .Examining a cross section of the distribution in Fig. 8b 7 at z = 5.1 × 10 7 m we find peaks at v z = ±1.3× 10 7 m s −1 , that is to say, in the same place as in Fig. 7h.The peak describes a semicircle in the z-v z plane for z > 5 × 10 7 m in Fig. 8b

Experiment 3
Experiment 3 started from the same initial state as the other two experiments.It followed experiment 1, but differed from it in that the voltage was brought down all the way to zero before returning to its initial value of 3 kV.Figure 9 shows the phase space density f e,M (µ , v z ) in experiment 3 at z = 5.1 × 10 7 m, which is a position close to the double layer on its high potential side.This is the same position that is shown for experiments 1 and 2 in Figs. 3 and 7, respectively.The results of experiment 3 are similar to those of experiment 1. Figure 9 shows the initial state and the distributions at t = 4 and 8 s, which is enough to show the differences between the two experiments.
At t = 4 s the region containing trapped electrons had shrunk to a small patch near (µ , v z ) = (0, 0).That this region did not disappear completely is due to fluctuations that created a minimum in the plasma potential of V p = −58 V at z = 4.93 × 10 7 m at t = 4 s.The speed of the precipitating particles was also lower here than in experiment 1.At the end of the simulation run, when t = 8 s, electrons had become trapped in the previously empty region at low µ 8 × 10 −12 A m 2 .The phase space density was somewhat higher in this region of Fig. 9c than in the corresponding region of Fig. 3h for experiment 1.The phase space density in the part of the trapped region above µ = 8 × 10 −12 A m 2 was lower at the end of experiment 3 than at the end of experiment 1, due to the former being more thoroughly emp-tied when the voltage was reduced to zero rather than only to 1 kV.
Figure 10 shows partial distribution functions f e,M (z, v z ) for electrons of magnetospheric origin in experiment 3. The division at µ = 8.0×10 −12 A m 2 is the same as that in Fig. 4, and like in Fig. 9 we have only included t = 0, 4, and 8 s.This figure serves to emphasise the differences between experiments 1 and 3, and it confirms that the trapping of low µ electrons and the emptying of the high µ region on the high potential side of the double layer were both more efficacious in experiment 3 than in experiment 1.
In Fig. 10a 4 electrons moving in the positive z direction were reflected at z = 4.9 × 10 7 m, whereafter they returned toward lower z values.This happened due to the shape of the potential profile at t = 4 s, which is shown by the blue curve in Fig. 5a.When the voltage was lowered all the way to zero, the double layer shifted polarity so that it at t = 4 s was accelerating electrons away from Earth.

Comparison between the experiments
The change in voltage first imposed on the system was V = −2 kV in experiment 1, V = −3 kV in experiment 3, and V = +2 kV in experiment 2. Figure 5a shows the plasma potential as a function of z for the initial state and the mean of the plasma potential during 0.1 s after the target voltages of the respective experiments were reached.In order to make a quantitative assessment of how the plasma assumes the changed voltage we define a reference point at z = 4.94 × 10 7 m, which is where the initial plasma potential was 1500 V, and it is shown by the left vertical dashed line in Fig. 5a.The right vertical dashed line is located at z = 4.99 × 10 7 m, and we shall take the double layer voltage to be the potential drop between those two vertical dashed lines.The positions of these lines were chosen to accommodate the double layer in all three cases, and the distance between them is therefore larger than the width of the double layer, which was approximately 180 km in the simulation of the initial state.In Fig. 5b the part of V that is assumed at low altitude -defined as z > 4.94×10 7 m -is shown together with the high altitude part, that is to say, z < 4.94 × 10 7 m.The part of V assumed by the double layer is also shown in the figure, and for reference the circle shows the initial plasma potential at the reference point, which was 1500 V.
In Fig. 5c the values shown in Fig. 5b have been normalised to the total V for the respective experiment, and the initial plasma potential at the reference point, shown by the circle, has been normalised to the initial total voltage.
In experiment 1, 67 % of V was assumed below the reference point; in experiment 2 that value was 68 %; and in experiment 3, 80 % of V fell below the reference point.Furthermore, it is seen that the part of V which was below the reference point was almost completely assumed by the double layer, and the plasma potential changed only a few volts over the distance from the double layer to the ionosphere.Therefore, that part of the plasma potential curves in Fig. 5a is flat on the scale shown.
Figure 11 shows the plasma potential at the end of the simulations as functions of z for the three experiments together with the initial plasma potential.Panel (a) shows the full z range and panel (b) shows the region around the double layer.In experiment 2, where the voltage excursion was toward higher values, the double layer was at the same position in the final as in the initial state.In the two other experiments, where the voltage was decreased, the double layer moved toward higher altitudes at the end of the simulation runs.This hysteresis effect was largest in experiment 3, which had the largest | V | of all three experiments.The double layer moved 370 km in experiment 3.
In all three experiments fluctuations in the plasma potential are seen in Fig. 11 at altitudes just above the double layer.This region extends to higher altitudes in experiment 2 than in the other two.In the region where these waves appeared, holes in the electron phase space can also be seen in Fig. 8 for t ≥ 5 s.
Figure 12 shows the phase space density of magnetospheric electrons and ionospheric ions in experiment 2 in panels (a) and (b) respectively for time t = 5 s.Panel (c) shows the plasma potential at t = 4.9 s (red line) and t = 5 s (black line), and panel (d) shows the electric field at these same times.In panels (a) and (b) the distribution function has been integrated over the complete µ range, and all panels show the z range 4.5×10 7 ≤ z ≤ 4.9×10 7 .The ion beam in the region shown was formed through acceleration in the double layer when the total voltage was increased to 5 kV and it was therefore faster here than at lower z values.
Each electron hole in Fig. 12a corresponds to a perturbation of the ion beam in Fig. 12b and a positive potential peak in Fig. 12c, showing that the hole was positively charged.The velocity with which the hole at z = 4.8 × 10 7 m moved was −7.5 × 10 5 m s −1 , determined from the zero crossings of the two curves in Fig. 12d.The hole was thus a manifestation of an ion beam mode -the ion beam velocity being in the same range, as seen in Fig. 12b.
To see whether the presence of the holes was caused by an ion-beam plasma instability, we have employed the method described by Gunell and Skiff (2001), where approximate Maxwellians are fitted to the beam and background ion distributions and the dispersion relation is computed using linear kinetic theory.The analysis shows that the beam mode is unstable and that the growth rate is approximately 0.06ω pi .This means that the typical growth time is many ion plasma periods.The simulations were run with r = 4.97 × 10 4 , and as the effective plasma frequencies scale as √ r the growth time is approximately 200 s at the location of the electron hole.

H. Gunell et al.: Auroral electrons trapped and lost
This is much longer than the duration of the experiments.Therefore, the excitation of the holes in the simulation could not be caused by the growth of this instability.In plasmas with two ion species, ion-ion instabilities would be active but on even slower time scales (Main et al., 2006).Instead the holes were created by the counterstreaming electron beams that are seen in Figs. 4, 8 and 10 that caused an instability on the much faster electron time scale.The positive net charge of the hole was the source of an electric field directed out of the hole, and that field then caused the perturbation of the ion beam, and the whole entity moved with the speed of the ion beam mode.The size of the hole at z = 4.8 × 10 7 m was on the order of 10λ e , where λ e = 47 km is the effective electron Debye length at the position of the hole.The effective Debye length is proportional to √ r ≈ 223, cf.Sect. 2. The real electron Debye length is therefore about 0.2 km and the expected size of such a hole in this region of space approximately 2 km.
Less conspicuous, but nevertheless similar, structures were seen in experiments 1 and 3.In the high µ range shown in Fig. 10b 8 , a hole was located at z = 4.8 × 10 7 m, and in the low µ range (Fig. 10a 8 ) there was a population of trapped electrons at the same z coordinate and low |v z |.A perturbation of the precipitating electrons also appeared at the same location.This was caused by the local potential peak at z = 4.8 × 10 7 m shown by the blue curve in Fig. 11b.

Conclusions and discussion
The auroral current circuit can be viewed as being composed of four elements: a generator in the equatorial magnetosphere and a load in the ionosphere that are coupled together by the upward and downward current regions.We model only the upward current region, which is where the electrons that cause the optical emissions are accelerated.This circuit element is affected by the other parts of the circuit.In order to examine what happens to the plasma on an auroral field line in response to fluctuations in the generator or the ionosphere, we have performed numerical experiments by changing the acceleration voltage during the run of a Vlasov simulation.We started from an initial state where a simulation with parameters of a typical auroral arc had been allowed to run until it reached a steady state.
We find that the electron distribution functions on the high potential side of the double layer, in the region between it and the ionosphere, are affected differently depending on whether the imposed voltage excursion is toward higher or lower voltages.In the case of a voltage excursion toward higher voltages -that is to say experiment 2 in this paper -the final distribution function f e,M (µ , v z ) was quite similar to the initial, although some structuring took place at intermediate values of µ .In the opposite situation, when the acceleration voltage was first lowered and then raised to its initial value, the distribution function changed significantly.
Electrons that had been trapped between the double layer and the magnetic mirror at high µ values were lost, and new electrons were trapped at low values of µ .
The way a particular voltage is reached has a marked effect on the shape of the distribution function.It should therefore be possible to say something about the recent history of the acceleration voltage by measuring the electron distribution function.A dominance of electrons at low µ is an indication of the voltage previously having been lower, and an abundance of high µ electrons indicates that the present voltage was reached from higher values.
The boundary between the high and low µ regions that behave differently is at µ = 8×10 −12 Am 2 .This is at the intersection between the two approximate parabolas determining whether the electrons can reach the ionosphere, given by Eq. ( 5), and whether they are able to reach the equator, which is determined by Eq. ( 6).When the total voltage is high, the electrons that start at the equator, and at z = z 1 = 5.1×10 7 m are in the region labelled II in Fig. 1g, are very close to the boundary that is given by Eq. ( 6), and a small perturbation may move them into region III where they become trapped.For electrons with lower µ values -that is to say, region I in Fig. 1g -the precipitating electrons are separated in velocity from the boundary given by Eq. ( 5), and here larger fluctuations are needed to move them into region III.The direction required to move the low µ electrons from region I to region III is toward lower v z , which corresponds to a lowering of the acceleration voltage, and this is consistent with trapping at low µ values occurring in experiments 1 and 3, where the voltage was lowered.
We found that when the voltage was changed, most of the change in voltage, V , was assumed by the double layer.A similar result was obtained by Forsyth et al. (2012), who used the Cluster spacecraft to show that most of V was taken up by the potential profile below the spacecraft trajectory.In our simulation the potential drop below the double layer altitude was no more than a few volts.To confirm this in a satellite observation one would need two spacecraft, one on either side of the double layer.If the total voltage decreases significantly and quickly, the polarity of the double layer is reversed for a period of time.In our experiment 3, where the voltage was lowered all the way to zero, the polarity of the double layer remained reversed until the voltage was increased again.Examining Fig. 11, we see that there is hysteresis in the double layer position for the experiments where voltage was decreased before it was increased.In both experiments the double layer moved toward higher altitudes.
In the simulations reported here, we changed the voltage in 1 s in experiments 1 and 2, and in 1.5 s in experiment 3 (see Fig. 2), and the complete simulations covered 7 and 8 s for these experiments respectively.The fluctuations observed by the Viking spacecraft were on the order of 1 Hz (André and Eliasson, 1992;Eliasson et al., 1996).The Cluster spacecraft, on the other hand, observed variations on a time scale of one or a few minutes (Marklund et al., 2001(Marklund et al., , 2011;;Sadeghi et al., 2011;Forsyth et al., 2012).The frequency of the fluctuations in Fig. 6 is about half a Hertz, and the time scale of the imposed voltage changes is a few seconds (Fig. 2).Thus, our simulations are closer to the Viking than the Cluster observations.However, the Cluster observations were limited by the spacecraft separation, and it is possible that the changes actually happened faster than it may appear.The time scale of the voltage changes that we imposed on the simulated system are similar to typical electron transit times reported in Sect.3.This situation is favourable for excitation of the kind of oscillations that are seen in Fig. 6.If the voltage rate of change had been slower, these oscillations likely would have been lower in amplitude and the reversal of the double layer electric field in experiment 1 less pronounced.However, if the voltage is brought down low enough, such as in experiment 3, the double layer remains reversed until the voltage is raised again or until the ions have had time to move out of the system.The latter alternative would take longer than the one minute time scale that the Cluster spacecraft observations set as an upper limit.
In experiment 2, where the voltage was raised to 5 kV before it was lowered again, Fig. 11 shows that a second double layer had developed at the top of the region where waves and phase space holes are present.Alm et al. (2014) found, based on Cluster satellite measurements, that the auroral acceleration cavity extends to high altitudes, and that 99 % of the potential drop is below these altitudes.Our simulations agree with that observation, the densities shown in Fig. 11c being more or less flat above the altitude of the double layer.That we find a cavity at high altitudes is no surprise, as we have used a cavity value for the density as a boundary condition on the magnetospheric side of the system.What the comparison shows is that such an assumption leads to solutions that are supported by space observations.Alm et al. ( 2014) furthermore concluded that "A midcavity double layer appears an unlikely explanation for the parallel potential structure".In light of the simulations conducted here, it seems possible that a midcavity double layer could exist at least for a limited period as a result of a varying voltage.However, this double layer only carries a small part of the total voltage.Ergun et al. (2004) observed mid cavity double layers that held a small fraction of the potential drop in the presence of wave turbulence.It is thus likely that midcavity double layers are dynamic in nature.
The width of the double layer reported in Sect.4.4 was 180 km.In this work we have used an artificial r = 4.97 × 10 4 for computational reasons.The effective Debye length is proportional to √ r ; the size of the double layer is expected to scale accordingly, and in space the double layer width would be about 1 km.As long as the double layer width is much smaller than the typical length scales of the overall changes of the plasma properties, which are about 5000 km on the high potential side and even longer on the low potential side, the exact width is not important for the results of the simulation.
The same scaling with √ r applies to the size of the electron holes shown in Fig. 12.This size was estimated to be about ten effective Debye lengths, which means that such holes in the auroral acceleration region should be about 2 km wide.The associated electric field is inversely proportional to the size of the hole and would therefore be expected to be about 1 V m −1 .Pottelette and Treumann (2005) observed tripolar electric field structures, which they interpreted as coupled electron and ion holes in accordance with simulations by Goldman et al. (2003).In the case of our electron holes, the electric field is bipolar, and the holes are formed by an instability of the counterstreaming electron beams that exist for a short period of time on the low potential side of the double layer (Figs. 4,8 and 10) as a result of the variable voltage that was imposed on the system.Once the holes have formed, they cause a perturbation of the ion beam, and they continue to exist even in the absence of the original electron beam responsible for their creation, moving away from Earth with the velocity of an ion beam mode.
The holes that are seen on the low potential side of the double layer in this study are not expected to be a source of radiation of any significance, as they are slow moving holes carrying no oscillating current.The holes that have been suggested as a source of auroral kilometric radiation appear on the high potential side, and are caused by the interaction of the electrons that are accelerated in the double layer and the background plasma (Pottelette et al., 2014).The emission mechanism they suggested is an electron cyclotron maser mechanism, where the electron distribution function is governed by the hole and may be different from the much studied horseshoe distribution (Ergun et al., 2000a).Similar systems have been seen to emit radio waves in laboratory plasmas (Lindberg, 1993;Volwerk, 1993;Brenning et al., 2006), although the electron beam energies were far below those used in dedicated cyclotron maser experiments (McConville et al., 2008) and simulations (Speirs et al., 2014).The simulations reported here are electrostatic and can therefore not be used to study radiation processes, which would require electromagnetic models.Electrostatic waves on electron time scales do appear on the high potential side of the double layer (Gunell et al., 2013a).The changes in the distribution function that occur as a result of the fluctuations studied in the present paper affect the beam-plasma interaction behind the excitation of these waves, and while a detailed wave study is beyond the scope of this paper it should be a topic for future research.

Figure 1 .
Figure 1.Initial state of the system for all three numerical experiments.(a) Plasma potential.(b) Electric field.(c)-(f) Phase space density f (z, v z ) for (c) magnetospheric electrons; (d) magnetospheric ions; (e) ionospheric electrons; and (f) ionospheric ions.(g)Phase space density f e,M (µ , v z ) for the electrons from the magnetosphere at z = 5.10 × 10 7 m, which is close to the double layer on its high potential side.The phase space regions separated by the white lines are I: precipitating electrons; Ib: up-going electrons that will reflect and then precipitate; II: electrons which can reach the equator; and III: trapped electrons.The unit for f e,M (µ , v z ) is m −6 A −1 s.In panels (c)-(f), the colour scales have been normalised so that integrals over all v z give n s /B.The unit for f (z, v z ) is m −4 T −1 s.

Figure 2 .
Figure 2. The total voltage imposed on the simulated flux tube as a function of time for the three different numerical experiments.

Figure 3 .
Figure 3. Phase space density f e,M (µ , v z ), in units of m −6 A −1 s, in experiment number 1 at z = 5.1 × 10 7 m for different times (on per second) are shown in the left-hand column.The time is indicated on each panel, and to the right that time is marked on a diagram of the total voltage applied over the system.

Figure 4 .
Figure 4. Partial phase space density f e,M (z, v z ) of magnetospheric electrons for experiment number 1. (a) Phase space density integrated over 0 ≤ µ < 8.0×10 −12 A m 2 .(b) Phase space density integrated over 8.0×10 −12 A m 2 ≤ µ < max(µ ) = 7.0×10 −9 A m 2 .The colour scales have been normalised so that integrals over all v z are equal to n p /B, where n p denotes the partial density of magnetospheric electrons for the µ range in question, and the unit for f e,M (z, v z ) is m −4 T −1 s.(c) Voltage as a function of time with the voltage for each row indicated by "×".The subscript of the panel label indicates the number of seconds for which the distribution function is shown in that panel.

Figure 5 .
Figure5.Distribution of V .(a) Plasma potential as a function of z for the initial state (green) and the mean of the plasma potential during 0.1 s after t = 1 s for experiments 1 (red) and 2 (blue) and t = 1.5 s for experiment 3 (blue), i.e. when the respective target voltages were first reached.The left vertical dashed line is located at z = 4.94 × 10 7 m where the initial plasma potential was 1500 V.The right vertical dashed line is located at z = 4.99×10 7 m.(b) The part of V below (red "+") and above (black "×") z = 4.94×10 7 m and the part across the double layer (blue "×"), all plotted versus the total V .The blue circle shows the initial potential at z = 4.94 × 10 7 m.(c) The same quantities as in (b), but with the vertical axis normalised to the total V .

Figure 6 .
Figure 6.Illustration of the effect of a time dependent potential profile on the particle distributions.(a) Phase space density f e,M (µ , v z ), in units of m −6 A −1 s, in experiment 1 at z = 5.1 × 10 7 m and t = 1.5 s.(b) Partial phase space density f e,M (z, v z ), in units of m −4 T −1 s, of magnetospheric electrons for experiment 1, integrated over 0 ≤ µ < 8.0 × 10 −12 A m 2 .(c) Plasma potential as a function of time for z = 4.92 × 10 7 m (black), z = 4.975 × 10 7 m (red), and at the ionosphere, i.e. z = 5.5 × 10 7 m, (blue).

Figure 7 .
Figure 7. Phase space density f e,M (µ , v z ), in units of m −6 A −1 s, in experiment number 2 at z = 5.1 × 10 7 m for different times (on per second) are shown in the left-hand column.The time is indicated on each panel, and to the right that time is marked on a diagram of the total voltage applied over the system.
7 .In fact, a semicircle of varying thickness is present in all the Fig. 8b panels.The final thickness and radius with the peaks at v z = ±1.3× 10 7 m s −1 for z = 5.1 × 10 7 m are a result of the outer -high energy -part of the semicircle, being released through the double layer toward lower z values after the voltage decrease.The thus released population is visible in Fig. 8b 4 around z = 4.9 × 10 7 m and v z = −1 × 10 7 m s −1 .It arrives on the low potential side at lower |v z | due to deceleration in the double layer.

Figure 8 .
Figure 8. Partial phase space density f e,M (z, v z ) of magnetospheric electrons for experiment number 2. (a) Phase space density integrated over 0 ≤ µ < 8.0×10 −12 A m 2 .(b) Phase space density integrated over 8.0×10 −12 A m 2 ≤ µ < 1.4×10 −11 A m 2 .(c) Phase space density integrated over 1.4 × 10 −11 A m 2 ≤ µ < max(µ ) = 7.0 × 10 −9 A m 2 .The colour scales have been normalised so that integrals over all v z are equal to n p /B, where n p is the partial density in the µ range in question, and the unit for f e,M (z, v z ) is m −4 T −1 s.(d) Voltage as a function of time with the voltage for each row indicated by "×".The subscript of the panel label indicates the number of seconds for which the distribution function is shown in that panel.

Figure 9 .
Figure 9. Phase space density f e,M (µ , v z ), in units of m −6 A −1 s, in experiment number 3 at z = 5.1 × 10 7 m for t = 0, 4, and 8 s are shown in the left-hand column.The time is indicated on each panel, and to the right that time is marked on a diagram of the total voltage applied over the system.

Figure 10 .Figure 11 .
Figure 10.Partial phase space density f e,M (z, v z ) of magnetospheric electrons in experiment number 3 for t = 0, 4, and 8 s. (a 0,4,8 ) Phase space density integrated over 0 ≤ µ < 8.0 × 10 −12 A m 2 .(b 0,4,8 ) Phase space density integrated over 8.0 × 10 −12 A m 2 ≤ µ < max(µ ) = 7.0×10 −9 A m 2 .The colour scales have been normalised so that integrals over all v z are equal to n p /B, where n p denotes the partial density of magnetospheric electrons for the µ range in question, and the unit for f e,M (z, v z ) is m −4 T −1 s. (c 0,4,8 ) Voltage as a function of time with the voltage for each row indicated by "×".The subscript of the panel label indicates the number of seconds for which the distribution function is shown in that panel.

Figure 12 .
Figure 12.Waves that are excited above the double layer in experiment 2. (a) Phase space density f e,M (z, v z ), in units of m −4 T −1 s, of magnetospheric electrons at t = 5 s.(b) Phase space density f i,I (z, v z ), in units of m −4 T −1 s, of ionospheric ions at t = 5 s.(c) Plasma potential.(d) Electric field.