The evolution of electron overdensities in magnetic fields

When a neutral gas impinges on a stationary magnetized plasma an enhancement in the ionization rate occurs when the neutrals exceed a threshold velocity. This is commonly known as the critical ionization velocity effect. This process has two distinct timescales: an ion–neutral collision time and electron acceleration time. We investigate the energization of an ensemble of electrons by their self-electric field in an applied magnetic field. The evolution of the electrons is simulated under different magnetic field and density conditions. It is found that electrons can be accelerated to speeds capable of electron impact ionization for certain conditions. In the magnetically dominated case the energy distribution of the excited electrons shows that typically 1% of the electron population can exceed the initial electrostatic potential associated with the unbalanced ensemble of electrons.


Introduction
In 1942 Alfvén attempted to explain the positions of the planets in the solar system [1]. He proposed that a neutral gas impinging on a magnetized plasma would cause significant amounts of ionization when the flow exceeds a threshold velocity, later known as the critical ionization velocity (CIV). The threshold velocity for the ionization enhancement is given by where eφ iz is the ionization potential of the neutral species in electron volts and m n is the mass of the neutral particle. Since the threshold velocity depends on the mass, species falling towards the Sun ionized at different Solar radii, seeding the planets. The CIV concept was later confirmed experimentally by Fahleson [2]. Brenning [3] and Lai [4] present extensive reviews of the literature. The CIV phenomenon has been implicated in many astrophysical contexts: solar abundances [5], interstellar filaments [6], plasma propulsion [7], comets [8], the Io-Jupiter system [9] and in its first instance, by Alfvén, to explain the seeding of planets in the solar system. Laboratory experiments have easily re-created and verified the CIV interaction [2,10]. Experiments carried out in the ionosphere have met with mixed results [11].
The idea is simple: rearranging (1) shows the neutral has kinetic energy equal to its ionization energy threshold. This energy, however, is insufficient for direct electron impact ionization since ionization occurs in the centre of momentum frame. Alfvén explained the CIV phenomenon in the following way [12]: neutrals, with kinetic energy eφ iz , streaming across a stationary magnetized plasma, perpendicular to the magnetic field, collide with ions resulting in momentum transfer or charge exchange. The collisions cause a region of unbalanced charge which is not neutralized immediately by the resdistribution of the electron population because the magnetic field restricts their transport. Regions of unbalanced negative charge can continue to grow, through the neutralization (effectively, the removal) of positive charge, until the potential reaches φ iz : too great for the ions to escape since the maximum speed the ions can gain is v c . The electrons are accelerated by this potential to speeds capable of ionization. Figure 1 shows a cartoon of the CIV effect: the impinging neutrals transfer momentum to the ions, the new trajectories of the ions move them away from the electrons, whose transport is much more inhibited by the magnetic field, the electrons are heated to higher energies by their self-electric field.
More detailed models have been proposed, most depending on the formation of unbalanced charge densities, but details of the structure and electron energization vary. The lower hybrid instability, also called the modified two-stream instability [13] is most commonly used to explain the electron heating. In general any process that creates a negative charge imbalance which heats the electrons to ionizing energies can be described as the CIV effect. There is no consensus on a single process responsible for the CIV interaction, several processes may be implicated for different conditions.
There have been several simulation studies of the CIV effect [14]- [16]. These simulations have used particle-in-cell [17] techniques in one and two spatial dimensions to verify the CIV effect. The ionization enhancement has been recreated in all cases and shown to depend on the type and functional form of the considered cross-sections [15,16].
A comprehensive investigation of the CIV effect, as a lower-hybrid instability, using simulation was carried out by Machida and Goertz [14]. A mostly neutral beam impinges on a stationary magnetized plasma. There is some pre-determined level of ionization in the beam resulting in beam ions. The density of the beam, stationary and thermal ions (ions resulting from elastic collisions with neutrals or diffused beam ions) are followed, as is the electron density. A uniform beam of neutrals is assumed and so no neutral properties are calculated 4 in the simulation. Three components are used to calculate the electron temperature evolution: resistive heating, due to elastic collisions with neutrals; collective heating, due to plasma waves; and collisional loses resulting from various excitations and ionization collisions.
Simulation studies thus far have shown the validity of the CIV mechanism but there have been no attempts to describe the dependence of the energy enhancement on the density and magnetic field strength. We aim to investigate the distribution of electron energies resulting from the acceleration under the self-field for different parameters. From this it may be possible not only to determine the optimum conditions for ionization but also for the creation of specific excited states.
In section 2, we describe the framework and implementation of a model to investigate the electron heating occurring during the CIV effect. Results from the simulation are shown in section 3 and discussed in section 4.

Model of electron energization
The CIV interaction can be broken into two problems on separate timescales: the ion-neutral timescale and the electron-electron timescale. On the ion-neutral timescale pockets of unbalanced negative charge are produced through charge exchange and momentum-transfer collisions between the ion and neutral species. The electrons in the pockets of unbalanced negative charge are accelerated by their self-field on the electron-electron timescale. The electron transport is restricted in the directions perpendicular to the magnetic field. This restricted transport increases the time the electrons take to neutralize the charge imbalance. The simulation of the electron-electron stage will be discussed here.
This stage is simulated by considering an ensemble of electrons in a magnetic field. Only the unbalanced electrons are simulated; the quasi-neutral background is assumed to be present but is not explicitly included in the calculation. There are two reasons for this choice: firstly, the inclusion of the background plasma would add to the computational load; secondly, including the background plasma would require the size of the electron pocket relative to the plasma density to be defined, adding an extra parameter. If the background plasma is strictly quasineutral then the contribution to the potential is zero. An initial population of electrons is allowed to evolve using single particle equations. Electron-electron interactions between the overdense population are of course encapsulated in the evolution of the common potential; background electron (and ion) collisions are neglected, though, as is any consequent evolution of the magnetic field.

Model equations
The equations required to simulate the evolution of the electrons are These equations are the momentum equation, Poisson's equation and the electric field, respectively. To carry out the numerical simulation the governing equations are nondimensionalized. The variables are mapped as follows: where the zero subscript denotes the scale factor and the non-dimensional variables are denoted by theˆnotation. The momentum equation (2) becomes Poisson's equation (3) becomeŝ Equation (4) becomeŝ We also define the non-dimensional kinetic energyκ, The non-dimensional parameters are In an effort to reduce the set of non-dimensional parameters such that the number of independent parameters is minimized we set: The variable α is a scaling factor, used to ensure that the values of the non-dimensional electric field (Ê) calculated in the simulation are similar in magnitude to the magnetic field strength (B) so that the values of P E and P B control the relative strengths of the electric and magnetic field components of the Lorentz force. Thus, Simplifying the non-dimensional parameters leaves a single free parameter P E . It may seem that α is also a free parameter, however, it is purely a scaling factor for the non-dimensional values and does not affect the evolution of the electrons. The factor α is useful as a technical aid to the numerical scheme. The parameter P E can be written as the ratio of the electron plasma frequency of the unbalanced electron population (ω p ) and the electron cyclotron frequency (ω c ),

Particle motion
We can see that the highest energy particles will be those transported parallel to the magnetic field direction. This is intuitive as the transport of electrons perpendicular to the magnetic field is inhibited. We can see this is mathematically justified by taking the dot product of the velocity and the acceleration (12) v · dv dτ = P Ev ·Ê + P Bv · (v ×B) This expression shows that the energy of a particle only changes whenv ·Ê = 0. A semi-infinite cylinder of charge is a good example containing bothv ·Ê = 0 and v ·Ê = 0 conditions. Figure 2 shows a diagram of a semi-infinite cylinder of charge annotated with electric field lines. Along the body of the cylinder the electric field is entirely perpendicular to the z-axis but at the end of the cylinder there is a fringing field and there are components of the electric field along each axis. We assume the magnetic field is orientated along the z-axis. The velocity of a particle gyro-orbiting the body of the cylinder will be perpendicular to the electric field and consequently will no gain energy sincev ·Ê = 0. A particle positioned at the end of the cylinder in the fringing field will always havev ·Ê = 0 even if it is gyro-orbiting the cylinder.

Numerical model
A three-dimensional mesh based Multigrid [18,19] solver is used to solve the finite difference representation of Poisson's equation (3) on a 129 3 grid. All calculations have been carried out with zero value Dirichlet boundaries, i.e. the potential at each boundary is fixed at zero. The plasma is assumed to be strictly quasi-neutral outside the computational domain, setting all boundaries to have a fixed potential of zero represents this assumption computationally.
When an electron reaches the boundary of the computational domain the particle is assumed to be absorbed by the quasi-neutral plasma and is removed from the simulation and plays no further part.
The particle trajectories are integrated using a fourth-order Runge-Kutta scheme for the ordinary differential equation of (12). The time step is chosen such that the motion of the particles is well resolved: in all cases the time step must resolve the Larmor orbits of the electrons but when the self-field dominates (P E 1) a finer time step may be required.
Collisions between electrons and other species have not been included. Clearly inelastic collisions will alter the energy gain of the electron population; this loss factor will be species dependent. By not specifically including such reactions, the calculations provide and initial guide as to where CIV processes may be significant. For example, ionization of xenon, from its ground state, requires electrons with energies above 12.12 eV, and so by examining the calculated electron energy distribution, conditions likely to produce enhanced ionization can be inferred. (Xenon in fact has a metastable state at 8.31 eV consequently a two stage ionization process would need only lower energy electrons.)

Initial conditions
The magnetic field is fixed to be uniform throughout the computational domain and is orientated solely along the y-axis. The unbalanced electron population consists of three million electrons, initially at rest, distributed normally in three dimensions around the point (N x /2, N y /2, N z /2). The normal distributions in each dimension have a variance of N x /3. The ideal nature of the initial electron distribution serves to illustrate the key concept of the role of the self-field versus the imposed magnetic field in the evolution of the electron energy distribution. The actual distribution of electrons in pockets of charge imbalance in a laboratory or astrophysical setting may be more complex, but the essential physics will remain the same. Table 1 gives a summary of the initial conditions used in the simulations.

Results
The evolution of the electron population was calculated for several values of the parameter P E , which describes the relative sizes of the electric and magnetic forces. Results from three values of P E are presented to highlight the behaviour of the electron overdensity in three different regimes, where the system is: magnetically dominated; electrically dominated; and when the forces are similar. Figures 3-5 show the evolution of the electron overdensity where P E = 10 −2 , 1, 10 2 , respectively. Snapshots were taken at various times through the simulation of the parameters: electron density, electric potential, mean energy and the distribution of kinetic  Figure 3. Snapshots of the evolution of the electron overdensity where P E = 10 −2 (magnetically dominated). The y-axis lies along the vertical and x-axis along the horizontal. Magnetic field is along the y-axis. An animation of the evolution can be seen in movie 1, available from stacks.iop.org/NJP/11/063001/mmedia. energy. Images of the spatial parameters: electron density, potential and mean energy are slices through the domain in the x-y-plane through the centre of the initial spatial distribution in the z-coordinate.
The time slices shown are at four regular intervals between the start and the point when the first electrons reach the boundary. Times are expressed in terms of electron cyclotron periods (τ c = 2π t 0 ). Non-dimensional units are used for all other quantities. Later, in section 3.5, the non-dimensional units are converted to physically meaningful units.

Magnetically dominated: P E = 10 −2
In the magnetically dominated regime the transport of the electrons perpendicular to the magnetic field is highly restricted and initial spherical electron distribution evolves towards a cylindrical shape. Figure 3 shows a series of snapshots from the simulation. 1. An animation of the evolution can be seen in movie 2, available from stacks.iop.org/NJP/11/063001/mmedia. Figure 4 shows that the evolution of the electron overdensity for P E = 1 appears to be very similar to the evolution when P E = 10 −2 . There are some differences most notably the timescale: when P E = 1 the electrons reach the boundary a factor of ten earlier. A peak in the kinetic energy distribution can also be seen.

Self-field dominated: P E = 10 2
For the case of P E = 10 2 (figure 5) the distribution of electrons evolves much quicker than the previous two cases. The electron density has dropped substantially after a quarter of a cyclotron period. This drop in density results in a flattening of the potential, in contrast to the previous two cases where the final potential was still relatively high.  Figure 5. Snapshots of the evolution of the electron overdensity where P E = 10 2 (self-field dominated). An animation of the evolution can be seen in movie 3, available from stacks.iop.org/NJP/11/063001/mmedia.

Comparing results
The ions that are displaced to form the initial charge imbalance will enter gyro-orbits and return approximately to their initial positions, neutralizing the charge imbalance after a gyroperiod. Once the charge imbalance is negated the acceleration of the electrons would end. In each case here the energization of the electrons occurs on a timescale much less than the ion gyroperiod.
The kinetic energy distributions for a range of values of P E are shown in figure 6. The distributions show the kinetic energy of the electrons as the first electrons reach the computational boundary. The maximum magnitude of the initial potential is denoted by a dashed, black line. In the simulations where P E > 1 the kinetic energies of the electrons do not exceed the maximum initial potential. Clearly the magnetic field plays a key role in the energy enhancement.
The fraction of the electron population exceeding the maximum initial potential is shown in figure 7. When the self-field dominates (P E > 1) there are many high-energy particles but this maximum energy does not exceed the initial potential. When the magnetic field dominates the distribution is peaked at a lower energy but has a high energy tail which exceeds the initial potential. The fraction of particles exceeding 50, 75 and 110% of the maximum initial potential are also plotted. To put this in context, if an overdensity with twice the potential of the ionization energy is created then for P E < 1 up to 13% of the electron population could exceed the ionization threshold.
We now define an efficiency factor, η: the ratio of maximum electron energy to the maximum of the modulus of the initial potential, (34) Figure 8 shows the variation of the efficiency for different values of P E . The efficiency factor is clearly dependent on the relative strength of the magnetic field; when the magnetic field dominates the maximum achievable energies are greater than the maximum initial potential.

Dimensional results
The results presented thus far have been in non-dimensional units. At this point we will convert to dimensionalized units. This is a simple procedure; the variables from the simulation are multiplied by their scaling factors as defined in section 2.  Figure 7. Fraction of electron population reaching energies in excess of the initial potential φ t=0 . Fraction of electrons exceeding 50, 75 and 110% of the initial potential are also shown. For P E = 1: 12% of the electrons exceed half of the initial potential, 3% exceed three-quarters of the initial potential and approximately 0.5% gain energies greater than or equal to the maximum initial potential. energy colour scale. Higher number densities result in higher maximum energies since the potential is related directly to the number density. For a given number density increasing the magnetic field strength gives a greater maximum energy.

Discussion
In this section we note some features of the electron overdensity evolution and in section 4.2 we attempt to explain these observations.

Characteristics
One of the most interesting features of the simulation results are the shapes of the electron kinetic energy distributions. The distributions all show similar morphologies: a peaked distribution with wings at higher or lower energies, or both (see figure 6). To investigate the origin of the distribution shape we separate it into high-and low-energy parts, using the peak of the distribution as the dividing point. Slices of the mean energy showing the high-and low-energy parts are shown in the bottom panels of figure 10, the kinetic energy distribution is shown in the top panel.
In regimes where the magnetic field plays a significant role, such that the distribution evolves over several cyclotron periods, the energy of the peak of the distribution can be seen to oscillate with a frequency approximately equal to the electron gyrofrequency. An example of this can be seen in the animation in movie 4, available from stacks.iop.org/NJP/11/063001/mmedia. The temporal evolution of the distribution of kinetic energies can be seen in figure 11.   Figures 3, 4 and 10 show that there is a preferred direction for the highest energies as expected in section 2.2. We can now bring the understanding from section 2.2 to the evolution of the electron distribution. Figure 12 is a sketch showing the electron gyromotion perpendicular to the magnetic field direction. The ensemble of electrons are initially in a spherical distribution. The electrons are accelerated by the self-field. The electrons located in the highest electric field positions gain a similar maximum energy and there is a sharp fall off in the kinetic energy distribution above the maximum energy. Electrons in other positions of the distribution gain lower amounts of energy filling out a low energy wing in the distribution. All electrons are travelling parallel to the electric field, i.e.v ·Ê = 0 and all electrons are gaining energy ( figure 12, point A).

Explanation
After a quarter of a cyclotron period (figure 12, point B) the electrons initially travelling at an angle to the magnetic field will be moving perpendicular to the electric field and the amount of energy being acquired by these electrons will be reduced sincev ·Ê ≈ 0 . This reduction in the gain of energy will effect a large proportion of the electrons including those represented by the peak in the kinetic energy distribution. The energy value of the peak will stop increasing. Some electrons will be moving parallel to the magnetic field; their transport will not be hindered by the magnetic field and they will continue to accelerate, creating a high-energy wing.
For the next half of the gyroperiod (figure 12, point C), those particles undergoing gyromotion will be travelling back across the magnetic field axis. This will most likely be travelling against the electric field and they will be decelerated. This reduction in acceleration can be seen in the oscillation of the peak in the kinetic energy distribution.

Astrophysical and laboratory context
It is worth commenting on the contrasting astrophysical and laboratory context of Alfvén ionization. The CIV phenomenon has been easily reproduced in laboratory experiments but space experiments have mostly failed to show any ionization enhancement. Comparing the experimental parameters [4], shown in table 2, with the results in figure 9 may explain the experimental discrepancy. Note that our electron densities refer to the overdense population only (that is, excluding any background plasma), and so the comparison between our results and the absolute electron density values given in table 2 may not be perfect, but it is strongly suggestive of the underlying reason for inconclusive astrophysical experimental results.
The conditions for the space experiments lie approximately in the bottom left of figure 9. Here the energy enhancement produced exclusively by acceleration under the electron self-field is less than 0.1 eV, substantially less than any ionization threshold. Of course, in our analyses we have not taken into account the streaming instabilities of the excess electron population and any background plasma; such instabilities are implicated in the generation of lower hybrid waves, considered to be a possible mechanism for ionization in such flows. However, such waves have proved elusive in observations [4,20], leaving the unbalanced self-fields caused by charge separation or ion removal as plausible seats of acceleration, and so the analysis in this paper is still relevant. In contrast, homopolar devices operate in the density and magnetic field strength regime located near the top right of figure 9. In this region, we expect electron energies to exceed ionization thresholds.
Hence one possible reason for the occurrence of the CIV effect in the laboratory but not in the space experiments is the charge density: without a sufficiently high plasma density, pockets of unbalanced charge with potential equal to the ionization potential cannot be created.

Conclusions
We have demonstrated, through numerical simulation, that electrons can be accelerated to energies capable of impact ionization by the self-field of a charge imbalance. The physical processes involved in the CIV interaction can be grouped into two timescales: an ion-neutral stage and an electron energization stage. This paper focuses on the short timescale energization part. The ion-neutral stage is responsible for the formation of regions in plasma where ions have been displaced leaving a negative charge imbalance. Calculating the evolution of these pockets of unbalanced electrons gives an estimation of the energy that can be obtained. Figure 9 shows that for certain charge imbalance densities and magnetic field strength values that ionization energies (typically φ iz > 10 eV) can be exceeded.
A single pocket of unbalanced charge consisting of three million electrons is simulated and the trajectory of each electron is calculated. The distribution of energies is important for determining the enhancement of the ionization rate. Figure 7 shows that when the magnetic field dominates the self field of the overdensity then approximately 1% of the electrons exceed the potential of the pocket. Under similar conditions more than 10% of the electron population exceed half the maximum initial potential.
The importance of the magnetic field was investigated by varying the parameter P E . The distribution of kinetic energies was found to be different in magnetically and electrostatically dominated regimes. Figure 6 shows the distributions for various values of P E . When the selffield dominates, the distribution is peaked and drops sharply towards high energy, consequently there is a large fraction of electrons in the high-energy part of the distribution but the energies do not exceed the maximum initial potential. The magnetically dominated cases (P E < 1) have wide high-energy tails reaching energies greater than the maximum initial energy. In this regime the transport of electrons perpendicular to the magnetic field is restricted; this means that the density structure evolves on a slower timescale and some electrons can gain energies in excess of the initial electrostatic potential.