Monte Carlo model for analysis of thermal runaway electrons in streamer tips in transient luminous events and streamer zones of lightning leaders

1.1. Streamers in Transient Luminous Events and Lightning Leaders

[2] Transient luminous events (TLEs) are large-scale optical events in the Earth's atmosphere that are directly related to electrical activity in underlying thunderstorms. Several types of TLEs are known: relatively slow-moving fountains of blue light, known as “blue jets,” that emanate from the top of thunderclouds up to an altitude of 40 km [e.g., Wescott et al., 1995, 2001; Lyons et al., 2003]; “sprites” that develop at the base of the ionosphere and move rapidly downward at speeds up to 10,000 km/s [e.g., Sentman et al., 1995; Lyons, 1996; Stanley et al., 1999]; “elves,” which are lightning-induced flashes that can spread over 300 km laterally [e.g., Fukunishi et al., 1996; Inan et al., 1997], and recently observed “gigantic jets,” which propagate upward, connecting thundercloud tops with the lower ionosphere [e.g., Pasko et al., 2002; Su et al., 2003].

[3] Since their discovery, there have been numerous imaging campaigns in an effort to better understand the physical phenomena behind these events. Recent imaging campaigns of sprites [Gerken and Inan, 2005; Marshall and Inan, 2005, and references therein], blue jets [Wescott et al., 2001], and gigantic jets [Pasko et al., 2002] have revealed a wide variety of fine filamentary structures in these events, which have been interpreted as streamers. Streamers are narrow filamentary plasmas, which are driven by highly nonlinear space charge waves [e.g., Raizer, 1991, p. 327]. Streamers can exhibit both positive and negative polarities, which is simply defined by the sign of the charge existing in the streamer head. Negative streamers generally propagate in the same direction as the electron drift, whereas positive streamers propagate opposing the electron drift. Negative streamers do not require ambient seed electrons to propagate since electron avalanches originating from the streamer head propagate in the same direction as the streamer [e.g., Vitello et al., 1994; Rocco et al., 2002]. Positive streamers, however, must obtain seed electrons from photoionization to sustain their propagation [e.g., Dhali and Williams, 1987; Raizer, 1991, p. 335].

[4] Streamers also serve as precursors to a more complicated leader phenomenon, which involves significant heating and thermal ionization of the ambient gas and represents a well known initiation mechanism of breakdown in long gaps [Raizer, 1991, p. 363]. Leaders are thin, highly ionized, highly conductive channels which grow along a path prepared by preceding streamers [Raizer, 1991, p. 364]. The head of the highly ionized and conducting leader channel is normally preceded by a streamer zone looking as a diverging column of diffuse glow and filled with highly branched streamer coronas [e.g., Bazelyan and Raizer, 1998, p. 203, 253]. The leader process is also a well-documented means by which conventional lightning develops in thunderstorms [Uman, 2001, p. 82], suggesting the presence of numerous streamers with every lightning discharge.

[5] It has been recently demonstrated that negative streamers developing in high ambient fields can reach an unstable “ideal conductivity” state with approximately equipotential and weakly curved head [Arrayas et al., 2002; Rocco et al., 2002]. This new state exhibits a Laplacian instability which can lead to branching of the streamer [Arrayas et al., 2002; Rocco et al., 2002] and can be realized over a wide range of applied electric fields [Liu and Pasko, 2004]. Liu and Pasko [2004] also studied the effects of photoionization on the dynamics of streamers and determined that the acceleration and expansion of streamers results in a reduction of the preionization level ahead of the streamers. In order to compensate for this reduction in preionization, the magnitude of the electric field in the streamer tip can reach a value as large as 10E_k at the stage immediately preceding the branching of the streamer, where E_k is the conventional breakdown threshold field defined by the equality of the ionization and dissociative attachment coefficients in air [e.g., Raizer, 1991, p. 135]. Figure 1a shows a negative streamer propagating at an altitude of 70 km in a 1.5E_k ambient field as it reaches an unstable state just prior to branching. It can be seen in Figure 1b that an extremely high electric field exists in the streamer tip prior to branching. This high field ∼2 kV/m, which spans approximately 1 m (see Figure 1b), could possibly accelerate low-energy electrons (∼several eV) to very high energies ∼2 keV. As discussed in the next section, the acceleration of electrons in these highly overvolted streamer tips could contribute to the formation of high-energy electron fluxes needed to explain the recently observed X-ray [Moore et al., 2001; Dwyer et al., 2003, 2004a, 2004b, 2005] and gamma ray [Fishman et al., 1994; Smith et al., 2005] bursts associated with thunderstorm activity.

[6] The model streamer shown in Figure 1 was obtained using the numerical model described by Liu and Pasko [2004] and assuming that no preionization is produced ahead of the streamer due to photoionization effects. These conditions are expected to be close to those realized before streamer branching when the photoionization range becomes shorter than the radius of the expanding streamer (see sections 4.1 and 4.4 in the work of Liu and Pasko [2004] for additional details).

1.2. Runaway Electrons and Energetic Radiation

[7] Runaway electrons were discussed by Gurevich [1961] and were defined by Kunhardt et al. [1986], who stated “an electron is runaway if it does not circulate through all energy states available to it at a given E/N, but on average moves toward high-energy states.” The runaway phenomenon is a result of decreasing probability of electron interactions with atomic particles for electron energies in the range from ∼100 eV to ∼1 MeV [Gurevich, 1961]. This phenomenon can best be understood by considering the dynamic friction force of electrons in air as a function of electron energy (see Figure 2):

where the summation is performed over all inelastic collision processes of a given gas with partial density N_j of N₂, O₂, or Ar in air (in m⁻³) corresponding to a particular collision process defined by the cross section σ_j and energy loss δε_j. In plotting F_D(ɛ) (Figure 2), electron-neutral collision cross sections provided by Phelps (http://jilawww.colorado.edu/www/research/colldata.html) and mass radiative and collision stopping powers [International Commission on Radiation Units and Measurements, 1984] were used. Electron energy losses due to nonzero energies of secondary electrons emerging from ionizing collisions with N₂, O₂, and Ar were accounted for using the differential ionization cross sections provided by Opal et al. [1971] (see 7).

[8] Electrons under the influence of an electric field E experience a force F_E = −q_eE and an acceleration according to the Lorentz force law and Newton's second law, respectively, where q_e is the absolute value of electron charge. As the electron accelerates through a gas it experiences collisions with the neutral gas molecules and atoms, which give rise to the dynamic friction force F_D opposing the force applied by the electric field F_E. It can be seen in Figure 2 that the friction force F_D varies considerably with electron energy. For example, a maximum exists in F_D at ∼100 eV which is ∼10³ the value of F_D at 1 eV. Physically speaking, a 100 eV electron experiences many more collisions and loses much more energy per unit length of its trajectory than does a 1 eV electron.

[9] F_D has units of eV/cm and can be directly compared to the applied electric field to provide an intuitively simple insight into the expected motion of electrons at various energies. Figure 2 lists electric fields required to initiate various types of electrical breakdown in air (more details and related references may be found in the work of Pasko [2006]) and displays the respective force F_E, in units of eV/cm, they apply to electrons in relationship to the friction force F_D. Of particular interest to the theory of runaway electrons is the maximum in F_D at ∼100 eV and the corresponding electric field which is known as the thermal [Gurevich, 1961] runaway threshold (E_c). Electrons with energies ∼100 eV moving through air will experience many collisions with neutral particles, which give rise to a high value of F_D. If an electric field E < E_c is applied to the electrons, the force F_E will be less than the force F_D the electrons will experience from collisions; therefore the electrons will be maintained at energies <100 eV. However, if an electric field E > E_c is applied to the electrons, it can be seen from Figure 2 that F_E > F_D. The electrons will gain more energy from the electric field than they will lose to collisions. It then becomes possible for some of the electrons to be energized to energies >100 eV. Owing to the reduced probability of collisions of electrons with energies >100 eV, the electrons will continue to accelerate to very high runaway energies as long as the electric field is present. Electric fields above E_c are difficult to produce and maintain since the electron runaway is also accompanied by an avalanche multiplication of electrons and strong increase in plasma conductivity, which tends to reduce the applied field. Electric fields ∼10E_k around tips of propagating streamers (Figure 1) are one of the unique naturally occurring circumstances when such high fields can be dynamically produced and sustained for relatively extended periods of time.

[10] Also of interest is the minimum in F_D which occurs at ∼1 MeV. This is known as the relativistic [Gurevich et al., 1992; Roussel-Dupre et al., 1994] runaway threshold (E_t) and is the basis of the Relativistic Runaway Electron Avalanche (RREA) model proposed by Gurevich et al. [1992]. At electron energies around 1 MeV the probability of collisions with neutrals is greatly decreased and any electron with an initial energy in this region (e.g., cosmic ray secondaries with energies 0.1–1 MeV [e.g., Roussel-Dupre et al., 1994; Gurevich and Zybin, 2005]) will run away when an electric field >E_t is applied. According to the RREA model, as few as one energetic electron (∼1 MeV) can trigger an avalanche of runaway electrons, via ionization of air molecules and atoms, which will continue to grow as long as an electric field E > E_t is present.

[11] Additionally, at lower electric fields comparable to the conventional breakdown field E_k, electrons are expected to be held to energies <20 eV by collisional losses since F_D > F_E for energies >20 eV. At even lower field values (i.e., fractions of E_k) electrons will be trapped by the local maximum in F_D around 1–2 eV resulting from strong energy losses due to excitation of vibrational degrees of freedom of nitrogen and oxygen molecules. The discrete structure observed at energies <1 eV also arises from the excitation of rotational degrees of freedom of nitrogen molecules and the excitation of rotational and vibrational degrees of freedom of oxygen molecules.

[12] The production of runaway electrons in the Earth's atmosphere has recently been linked to X-ray and gamma ray bursts observed during lightning discharges [Moore et al., 2001; Dwyer et al., 2003, 2004a, 2004b, 2005]. In addition to these ground-based measurements, intense gamma ray flashes originating from the Earth's atmosphere above thunderstorms have also been observed by the Compton Gamma Ray Observatory (CGRO) and the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) satellites [Fishman et al., 1994; Smith et al., 2005]. While these observations strongly support the existence of extremely high-energy electrons during thunderstorm activity, the exact mechanism producing them remains under debate [e.g., Dwyer, 2005a]. Gurevich [1961] showed that in the presence of extremely strong electric fields, a large number of low-energy electrons can be directly accelerated over the peak of the friction force F_D and become thermal runaway electrons. This is a relatively straightforward approach to runaway development and is readily accepted. However, since the electric field strengths necessary to achieve thermal runaway (E ∼ 10E_k) and even conventional (E ∼ E_k) breakdown are not commonly observed on large spatial scales in thunderclouds [Marshall et al., 1995, 2005], many scientists abandoned thermal runaway breakdown as a source of runaway electrons during thunderstorms and adopted the newer theory of RREA. One of the goals of the present paper is to demonstrate that streamers may represent a realistic source of thermal runaway electrons and discuss circumstances when these electrons can be accelerated to very high (>1 MeV) energies, thus providing an alternate source of seed electrons to the RREA model previously thought to require galactic cosmic rays.

1.3. Purpose of This Paper

[13] This paper presents the formulation of a Monte Carlo model, which is capable of describing electron dynamics in air including the electron thermal runaway phenomena under influence of an external electric field of arbitrary strength. The model is similar in technical details to the model previously developed for N₂ by Tzeng and Kunhardt [1986] and incorporates the following features: (1) the “null” collision method to determine time between collisions [Lin and Bardsley, 1977]; (2) the remapping of the electron assembly to improve statistics for the high-energy tail of the electron distribution [Kunhardt and Tzeng, 1986b]; (3) the differential ionization [Opal et al., 1971] and scattering cross sections for realistic description of energy spectrum of secondary electrons and the forward scattering properties of electrons at high energies. Results for zero-dimensional modeling of the electron distribution under the influence of a uniform electric field are first presented and compared with existing data. At high electric fields the model is validated by comparisons with studies conducted for N₂ by Tzeng and Kunhardt [1986] and more recently by Bakhov et al. [2000]. At low electric fields, model results are compared to available data from swarm experiments in air [Davies, 1983], numerical solutions of the Boltzmann equation based on the two-term spherical harmonic expansion of the electron distribution function [Morgan and Penetrante, 1990], and analytical models proposed by Aleksandrov et al. [1995] and Morrow and Lowke [1997]. The new Monte Carlo model is then applied to a one-dimensional case corresponding to a negative streamer tip immediately preceding branching. The model results demonstrate that the electric fields in streamer tips are strong enough to accelerate low-energy electrons to several keV, initiating thermal runaway in relatively low ambient fields. Streamers have been documented in transient luminous events above thunderclouds and in streamer zones of conventional lightning leaders and may provide a robust source of runaway electrons contributing to the production of recently observed X-ray and gamma ray bursts.

2.1. Collision and Scattering Cross Sections

[14] Essential to both the numerical Boltzmann equation solution and the Monte Carlo method describing behavior of electrons in weakly ionized air is the knowledge of electron-atom and electron-molecule collision cross sections for each gas species. Generally, atoms and molecules in a weakly ionized gas are assumed to be heavy and remain stationary; therefore the collision cross section becomes simply a function of the electron energy.

[15] For the Monte Carlo model presented in this paper several different cross section sets were used to determine an electron's motion through air. A total of 46 collision processes were considered corresponding to the N₂, O₂, and Ar gases included in the air mixture (see Tables 1, 2, and 3). Each gas has an associated elastic cross section accounting for elastic collisions and a set of inelastic cross sections (24 for N₂, 17 for O₂ and 2 for Ar) accounting for each inelastic collision process (i.e., rotational, vibrational, and electronic excitations, ionization, and attachment).

[16] The 43 inelastic cross sections σ_inel were obtained directly from the compilation of A.V. Phelps (http://jilawww.colorado.edu/www/research/colldata.html). The elastic collision cross section σ_el for each gas was then determined by subtracting the summation of inelastic cross sections from the total collision cross sections

where s represents a specific gas species and the summation is performed over all inelastic collision processes with j representing the jth process. The total collision cross sections were obtained from experimental data reported in literature as summarized in Table 4.

[17] The calculated elastic cross sections and the inelastic cross sections are the fundamental quantities used to determine an electron's interaction with the gas medium. From these collision cross sections the mean free path, mean time between collisions, and collision frequency of an electron in the gas medium can be calculated as

respectively, where N is the gas density, σ is the collision cross section, and v is the electron's velocity.

[18] In addition to the collision cross sections mentioned above, the differential scattering cross section for each gas must be known to determine the angular scattering of electrons after a collision. Experimental values for the elastic differential scattering cross sections were obtained from literature (see Table 5). The angular scattering of electrons in inelastic collisions was determined using the elastic differential cross sections listed in Table 5. This assumption is reasonable for scattering of electrons on the molecules with excitation of singlet states (predominantly forward) but may introduce error for the scattering involving excitation of the triplet states, which possess a stronger backscatter component [Tzeng and Kunhardt, 1986]. Figure 3 shows the elastic differential scattering cross section, for low-energy scattering obtained from sources listed in Table 5 for N₂, O₂, and Ar. It can be seen from Figure 3 that at electron energies <20 eV the angular scattering of electrons is generally isotropic but quickly becomes forward for electron energies >20 eV. Electrons with energies ≈0 eV were assumed to demonstrate isotropic scattering and the differential cross section values from 0 to 5 eV for N₂, 0 to 2 eV for O₂, and 0 to 3 eV for Ar were determined using linear interpolation.

[19] For use in the Monte Carlo calculations, differential cross section tables for each gas were generated using linear interpolation of the data from sources listed in Table 5. The tables are then related to a uniform random number R_χ from 0 and 1 as

where χ is the scattering angle of an electron. The scattering angle χ can then easily be found by performing a table lookup using a random number R_χ and the electron energy ɛ. Angular scattering of electrons with energies in the range 0–500 eV experiencing collisions with N₂, and in the range of 0–1000 eV colliding with O₂ and Ar is determined using lookup tables derived from published experimental data from sources given in Table 5.

2.2. High-Energy Electron Scattering

[20] For collisions involving electrons with energies greater than those tabulated from the experimental data (>500 eV for N₂ collisions, >1000 eV for O₂ and Ar collisions), differential scattering cross sections are calculated using several different approximations. It is the shortage of cross section data at high electron energies that remains the largest source of error in Monte Carlo simulations and there is no commonly accepted approximation which is universally used.

[21] For low-energy electrons, the dominant process is relatively short-range polarization scattering [Lieberman and Lichtenberg, 1994, p. 60]. However, the mean collision time (equation (4)) of high-energy electrons is small, not allowing atoms and molecules adequate time to polarize. Therefore high-energy electron scattering from neutral particles (i.e., electron-atom and electron-molecule collisions) resembles Coulomb-like collisions between charged particles (i.e., electron-electron, electron-ion, and ion-ion collisions) [Lieberman and Lichtenberg, 1994, p. 60].

[22] The differential cross section for Coulomb-like collisions can be analytically derived as

where q₁ and q₂ are charges of projectile and target particles, ɛ₀ is the permittivity of free space, ɛ is the center-of-mass energy, and b₀ is known as the classical distance of closest approach. Equation (7) is the well known Rutherford cross section for Coulomb scattering [Lieberman and Lichtenberg, 1994, p. 57]. In a general case, however, the Rutherford cross section cannot be directly used in Monte Carlo simulations because of the singularity as χ → 0. For this reason, a wide number of approximations based on screened Coulomb scattering have been implemented to determine the angular scattering of high-energy electrons in Monte Carlo and Boltzmann equation solutions.

[23] A. V. Phelps (http://jilawww.colorado.edu/www/research/colldata.html) presents an analytical differential scattering cross section approximation for elastic electron scattering from N₂ based on a screened-Coulomb type scattering

where here and in subsequent equations ɛ is in units of eV and β(ɛ) is an algebraic screening parameter derived to fit experimental angular distributions from Phelps and Pitchford [1985]. Substituting equation (9) into equation (6) and solving for the scattering angle χ results in

where R_χ is a uniform random number between 0 and 1.

[24] Surendra et al. [1990] proposed an analytical expression based on screened Coulomb scattering from Ar

such that the electron-neutral scattering at low energies is mainly isotropic and becomes increasingly anisotropic as the electron energy increases. Substituting equation (12) into equation (6) and solving for the scattering angle χ results in

[25] In his study of an electron avalanche development in neon, Shveigert [1990] used the differential scattering cross section for the scattering of fast electrons by shielded Coulomb potential of a nucleus as published by Kol'chuzhkin and Uchaikin [1978]:

where Z is the number of protons in the atom's nucleus and η is the shielding parameter formulated to fit the calculations of Thomas [1969] at ɛ = 100 eV. Scattering angles from this approximation are tabulated numerically using equation (6).

[26] A modified Rutherford cross section was also introduced for use in this paper of the form

where ɛ₁ is a shape parameter which is set to ɛ₁ = 4 eV to match experimental data of Kambara and Kuchitsu [1972] for electron scattering from N₂ at 500 eV. Substituting equation (16) into equation (6) and solving for the scattering angle χ results in

[27] The differential cross section used to approximate high-energy electron scattering can greatly impact the generation of runaway electrons and will be discussed further in 20, 22, and 23. Figure 4 shows the differential cross sections described by equations (9), (12), (14), and (16).

[28] It should be noted that while Figures 4a through 4d are plotted for the energy range 0–1000 eV, the cross sections are only used in the model for collisions with N₂ where ɛ > 500 eV and collisions with O₂ and Ar where ɛ > 1000 eV. Figure 5a shows a comparison of equations (9), (12), (14), and (16) with experimental data of Kambara and Kuchitsu [1972] for electron collisions with N₂ at ɛ = 500 eV and Figure 5b shows a similar comparison with experimental data of Iga et al. [1987] for electron collisions with O₂ at ɛ = 1000 eV. Upon first examination of Figures 4 and 5, the differences between the differential cross sections may appear to be small. However, the forward scattering properties of high-energy electrons are vital in the development of electron runaway and these small variations can drastically hinder or facilitate the runaway process (20, 22, and 23).

2.3. Differential Ionization Cross Section

[29] The behavior of the phase-space distribution of electrons is strongly influenced by ionization and the angle and energy of each of the two electrons emerging from the collision [Tzeng and Kunhardt, 1986]. Tzeng and Kunhardt [1986] placed special emphasis on the energy partitioning used in ionizing collisions and presented results for four separate cases:

[30] In case 1 the secondary electron is assigned zero energy, leaving the primary electron with the difference between initial energy and the ionization energy. In case 2 the energy of the secondary is determined from the differential cross section for ionization determined from experiments by Opal et al. [1971]. In case 3 the fraction of half the available energy given to the secondary is a random variable uniformly distributed in the interval [0,1]. In case 4 the primary and secondary electrons equally share the available energy.

[31] The differences in the energy distribution functions for each of these cases are significant, especially at high electron energies. The variations can best be explained by considering the energy of the incident electron after an ionizing collision

where ɛ_p is the energy of the incident electron before the collision, ɛ_s is the energy of the secondary electron generated by the collision, and ɛ_iz is the ionization potential. For cases 3 and 4 the high-energy incident electrons participating in the ionizing collisions lose a large fraction of their energy to the secondary electron (ɛ_s is large), thus severely reducing the number of electrons which can accelerate to runaway energies. In contrast, incident electrons in case 1 lose only energy equal to the ionization potential, ɛ_iz, (ɛ_s = 0), thus contributing to an overpopulation of electrons existing at very low energies due to the low-energy secondaries as well as an overpopulation of runaway energy electrons due to high-energy incident electrons losing only a small fraction of their energy to ionizing collisions. Among the cases presented, case 2, which makes use of the secondary energy distribution suggested by Opal et al. [1971], is the most realistic [Tzeng and Kunhardt, 1986] and is used for all simulations presented in this paper.

[32] Opal et al. [1971] measured a quantity proportional to the doubly differential cross section σ(ɛ_p, ɛ_s, χ). Integrating this cross section over the angle χ, results in

Assuming the ion to be massive and at rest, the kinetic energy imparted to the ion in the collision is negligible and the energies of the two departing electrons must sum to ɛ_p − ɛ_iz and be symmetrical about (ɛ_p − ɛ_iz)/2 [Opal et al., 1971]. The total ionization cross section can then be given by

From results of their experiments, Opal et al. [1971] determined the differential ionization cross section to be

where is a shape parameter adjusted to fit the ejected electron spectrum. Values of determined by Opal et al. [1971] are listed in Table 6 for N₂, O₂, and Ar gases.

Table 6. Ionization Energies _iz and Ejected Electron Spectrum Shape Parameter in eV [Opal et al., 1971]

Gas	ɛiz
N2	15.6	13.0
O2	12.2	17.4
Ar	15.7	10.0

[33] The differential ionization cross section of Opal et al. [1971] can be used to determine the energy ɛ_s of the secondary electron. Similarly to the determination of the scattering angle from the differential scattering cross section from equation (6), ɛ_s can be related to a uniform random number by performing the integration

Observing that the denominator of equation (22) is the total ionization cross section σ_i(ɛ_p) and substituting the differential ionization cross section of Opal et al. [1971](21), equation (22) can be rewritten as

Solving equation (23) for the secondary electron energy ɛ_s results in

[34] Also, as mentioned in 2, the differential ionization cross section of Opal et al. [1971] was used in the determination of the dynamic friction force in air (see Figure 2). Using the differential cross section defined by equation (21), the average energy of a secondary electron emerging from an ionizing collision can be found as

After obtaining the average secondary energy 〈ɛ_s(ɛ_p)〉, the friction force of ionizing collisions can be calculated as

where the index j accounts for differences in the ionization potential and the average secondary energy corresponding to different target species with density N_j.

2.4. Null Collision Method

[35] The Monte Carlo method works by individually tracking each electron in an assembly through a series of time steps until an equilibrium state is attained (e.g., electron mean energy remains constant). In a given time step, Δt, the electron may or may not experience a collision with the probability of a collision being [Dincer and Govinda Raju, 1983]

where τ is the mean time between collisions of an electron defined by equation (4). We note that the mean collision time in equation (4) depends on the electron's velocity. Therefore for a large number of electrons included in a simulation, there will be an equally large number of distinct time steps between collisions. The computation required to support each electron possessing its own mean collision time is a daunting task and can be overcome by adopting a null collision technique [Lin and Bardsley, 1977], which allows a constant mean collision time to be defined for all electrons in the system. The null collision approach first defines the total collision frequency as the sum of all elastic and inelastic collision frequencies

where

m is the mass of an electron, and N_j represents a partial density of target molecules or atoms corresponding to a particular collision process defined by the cross section σ_j. Having plotted ν_t(ɛ) versus electron energy (as schematically shown at the top of Figure 6), the maximum collision frequency ν_max can be found

The constant mean collision time τ_c can then be calculated by substituting ν_max into equation (5)

A constant time step Δt can then be defined as

where δ is an arbitrary number much less than 1. In all calculations presented in this paper, δ is assumed to be equal to 0.1 (see Figure 6).

[36] Now, having obtained a constant time step to be used throughout the simulation, the null collision technique can be viewed as a three step procedure using three uniform random numbers between 0 and 1 to determine if an electron experiences a collision in a time step, and if so, what type of collision it was. An outline of the null collision method is provided in Figure 6.

[37] 1. Substituting τ_c and Δt into equation (27) results in a constant probability of a collision in a time step

Therefore if a uniform random number R₁ from 0 to 1 is less than 0.1 for a given electron, the electron is said to have experienced a collision in the time step.

[38] 2. Having determined that an electron experienced a collision in step 1, we must now determine whether the collision is a null or real collision. The energy independent maximum collision frequency ν_max defined by equation (30) can be represented as the sum of the energy dependent null ν_null(ɛ) and real ν_t(ɛ) collision frequencies (see top of Figure 6)

The probability of a real collision, P_real, can then be defined as

and is a function of the electron energy ɛ. For each electron which was determined to have experienced a collision in step 1, a second random number R₂ with a uniform distribution between 0 and 1 is generated and compared with the probability P_real corresponding to the electron's energy. If R₂ ≤ P_real the collision is considered to be real, otherwise (R₂ > P_real) the collision is considered to be null and has no effect on the electron's properties (see Figure 6).

[39] 3. If the collision was determined to be real in step 2, then the next step is to determine what type of collision (i.e., momentum transfer, excitation, ionization) occurred. At a given electron energy ɛ_j, there is a collection of individual collision frequencies ν_j(ɛ), which sum to equal the total collision frequency ν_t(ɛ) as shown by equation (28). The probability of each collision process can then be calculated as

at a given electron energy ɛ such that

Each collision process can then be assigned a range of numbers existing between 0 and 1 by performing a cumulative summation of the individual probabilities P_j such that Range 1 = 0 to P₁, Range 2 = P₁ to P₁ + P₂, Range 3 = P₁ + P₂ to P₁ + P₂ + P₃, etc. A third uniform random number R₃ between 0 and 1 can then be generated and whichever process's range it falls within is determined to be the collision process which occurred (see Figure 6).

[40] To better illustrate the procedures of step 2 and step 3, consider a simple case of only argon gas. The collision frequencies of argon are shown in Figure 7a. Normalizing the collision frequencies according to equation (35) results in values between 0 and 1 as shown in Figure 7b. From Figure 7b it can be seen that for a 100 eV electron, the total normalized collision frequency is 0.0825 + 0.3092 + 0.4709 = 0.8626. A real collision is said to occur if random number R₂ is such that 0 < R₂ < 0.8626 and a null collision occurs if R₂ > 0.8626. If it is determined that a real collision has occurred, then the collision processes must be normalized by the value of the total collision frequency at 100 eV of 0.8626 as shown in Figure 7c. Now random number R₃ will be generated to determine which type of collision occurred. If 0 < R₃ < 0.5459, the collision is said to be elastic. If 0.5459 < R₃ < 0.6415 (where 0.5459 + 0.0956 = 0.6415), the collision is said to be excitation. Finally, if 0.6415 < R₃ < 1 (where 0.6415 + 0.3585 = 1), the collision is said to be ionization. After the type of collision is determined, the electron's properties (i.e., energy, velocity, direction) are modified accordingly as described in 9 below.

[41] After simulating each electron's interactions through a time step as described in steps 1, 2, and 3, the electron velocities are updated to reflect acceleration due to the applied electric field and diagnostic data is saved. If a specified number of electrons has been reached due to ionization processes, a particle remapping scheme (14) is applied to avoid undesirably long computation times. The procedure above is then repeated until an equilibrium state is reached or until a point in time, which is defined by an investigator.

2.5. Energy Loss and Scattering of Electrons

[42] There are four types of collision processes considered in this paper: elastic, excitation, ionization, and attachment. When one of these collision processes occur in a Monte Carlo simulation (see 8), the colliding electron's energy and trajectory following the collision must be determined. A summary of these collisions is shown in Figure 8.

[43] First, consider an electron characterized by its perpendicular v_⊥ and parallel v_∥ velocity components with respect to the applied electric field and angles θ_c and ϕ_c as shown in Figure 9a, where

Since only electron collisions with massive, stationary neutral atoms and molecules are considered, scattering events can be treated in center-of-mass coordinates [Lieberman and Lichtenberg, 1994, p. 51]. For simplicity, the initial angle ϕ_c is assumed to be zero and the electron velocity v is always in the (x, z) plane as shown in Figure 9b. Now envision the electron colliding with neutral particle resting in the line of v and scattering at an angle χ through a differential solid angle dΩ = sinχdχdϕ. To find the new trajectory of the electron after the collision, it is useful to define a new coordinate system (x′, y′, z′) by rotating the initial coordinate system (x, y, z) about the y-axis by θ_c as shown in Figure 9c. The electron scattering can now easily be treated in the new (x′, y′, z′) coordinates as shown in Figure 9d where v_old is the particle velocity before scattering, θ_c is the angle of the particle before scattering in (x, y, z) coordinates, v_new is the velocity after scattering, χ is the angle after scattering in (x′, y′, z′) coordinates, θ is the angle after scattering in (x, y, z) coordinates, and ϕ is considered to be random. The new trajectory of the electron in (x, y, z) coordinates can then be calculated as follows:

where cos θ_c and sin θ_c describe the electron's trajectory before the collision, cos θ and sin θ describe the trajectory after the collision, and R_ϕ is a uniform random number between 0 and 1.

2.5.1. Elastic Collisions

[44] In the case of an elastic collision, the only energy loss mechanism is the momentum transfer between the electron and the neutral particle. This energy loss is a function of the scattering angle of the electron following the collision, therefore the first step when an elastic collision occurs is to determine the scattering angle χ from equation (6). After obtaining χ, the new trajectory of the electron can be calculated using equations (41) and (42) and the energy loss due to momentum transfer can be determined as [e.g., Liu and Govinda Raju, 1992; Lieberman and Lichtenberg, 1994, p. 54]:

where M is the mass of the neutral particle.

2.5.2. Excitation Collisions

[45] Although the physics behind excitation collisions is complex, calculating the incident electron's properties after an excitation collision is trivial. The excitation energies and collision cross sections corresponding to many electron-neutral collisions have been well studied and experimental data is readily available in literature. Therefore the new energy of an incident electron after an excitation collision can be simply calculated as

where δɛ_j is the energy lost by the electron for the excitation process j, as schematically illustrated in Figure 8 for a case of excitation collisions. After obtaining the electron's new energy ɛ′_p, the scattering angle of the electron can be determined from equation (6) noting that ɛ is replaced by ɛ′_p. This approach is consistent with that used by Kunhardt and Tzeng [1986b]. The electron's new trajectory can then be calculated using equations (41) and (42).

2.6. Particle Remapping

[48] A limiting factor associated with Monte Carlo simulation is the long computation times and computer memory required to fully describe a physical system with millions of individual elements and a large number of time steps required to reach a converging solution. In the case of using the Monte Carlo technique to calculate the electron energy distribution in air, this problem arises due to increasing ionization rates at high electric fields leading to an enormous multiplication of particles to be tracked in the simulation. To resolve this issue, a so-called nonanalog Monte Carlo technique of statistical weighting similar to that of Kunhardt and Tzeng [1986b] is introduced.

[49] After a predetermined number of particles, N_t (usually 15,000), is reached in a simulation, a remapping of the electron assembly is performed to reduce the total number of particles to another predetermined number, T_t (usually 8000), with each new particle being assigned a new particle weight to reflect how many electrons the particle represents. Given the nature of the electron energy distribution function in air (i.e., the bulk of the electron distribution is maintained at lower energies with a small number of electrons constituting a high-energy tail), measures must be taken to maintain appropriate resolution in the high-energy tail of the distribution especially in the study of runaway phenomena. To do this, the initial particles N_t are first sorted in order of increasing energy. After sorting the particles, they are then partitioned into two groups, one group representing the low-energy body N₁ of the distribution and another representing the high-energy tail N₂ (see Figure 10). The first low-energy group of particles is remapped such that the weights of adjacent particles in energy space (W_j and W_j+1) are combined into one particle with a new particle weight W′_j such that

The new particle assumes the properties of the jth particle such that ɛ_new = ɛ_j, v_∥,new = v_∥,j, and v_⊥,new = v_⊥,j, which is an adequate assumption since the energy difference between two adjacent particles in energy space after the sorting will be small in the low-energy body of the electron energy distribution. This results in the total number of particles in the low-energy region being reduced from N₁ to T₁ (usually 7000) particles. The high-energy tail group is then remapped 1:1 to ensure enhanced resolution at high electron energies such that

for all particles in the tail and N₂ = T₂ (usually 1000). The procedure is repeated each time the particle count in a simulation exceeds N_t and is validated by comparing the electron energy distribution function from before and after the remapping events.

[50] For simulations when the spatial distribution of electrons must also be considered, as with the one-dimensional treatment of runaway electrons in a streamer tip (22), remapping the particle assembly becomes slightly more complicated. When remapping occurs, the particles are first sorted according to increasing position along the z-axis. The particles are then divided into equally spaced spatial bins (usually 10) and then sorted in order of increasing energy within each bin. The remapping of particles is then performed within each bin according to the same 1:2 for low-energy and 1:1 for high-energy particles as discussed in the previous paragraph (see Figure 11). The remapping scheme is validated by comparing the electron energy distribution function before and after a remapping event.

2.7. Model Initialization and Execution

[51] Figure 12 shows a flow chart representing the execution of the Monte Carlo model. First, the user must input basic simulation parameters such as the initial number of particles to be used (usually 8000), the applied electric field strength, the length of the simulation, the fractional composition of the gas mixture (i.e., 78.11% N₂, 20.91% O₂, and 0.98% Ar), and the number of particles at which the particle assembly will be remapped (usually 15,000, see 14).

[52] After the input parameters are entered, the model can then calculate the initial electron distribution and define necessary quantities to be used throughout the simulation such as the collision frequency and the size of the time step. The initial electron set is normally assigned a Maxwellian velocity distribution function [Chen, 1984, p. 226; Birdsall and Langdon, 1991, p. 390] corresponding to an initial electron temperature of T_e = 5800 K (i.e., 0.5 eV). To achieve a Maxwellian velocity distribution, each electron is assigned a v_x, v_y, and v_z velocity component. Using the error function, defined by

where in the context of the present problem 0 ≤ v ≤ 5 is a generic range normalized by thermal velocity v_th = , the normalized v_x component of an electron's velocity can be related to a uniform random number R_j between 0 and 1 from the relationship [e.g., Birdsall and Langdon, 1991, p. 390]

Similarly, the normalized v_y and v_z velocity components can also be found for each electron using equation (48). After assigning velocity components to half of the electron population, the remaining electron velocities may be found by mirroring the positive velocity components found from equation (48) to the corresponding negative values. A true three-dimensional (3-D) Maxwellian distribution is then arrived at by multiplying the normalized v_x, v_y, v_z components of each electron by the thermal velocity v_th. The 3-D distribution is then transformed to parallel and perpendicular velocity components for use in the Monte Carlo simulation as

assuming an electric field in the z-direction (see Figure 9).

[53] The total collision frequency ν_t(ɛ), the collision frequency for each process ν_j(ɛ), and the maximum collision frequency ν_max are calculated from the cross section data sets (5) prior to the simulation to reduce computation time. These values are then loaded into the simulation at runtime and used to calculate the time step (equation (32)) and to determine if and what type of collision an electron experienced in a time step as outlined in Figure 6. In addition, the differential scattering cross section data discussed in 5 is also tabulated prior to execution, allowing for convenient table lookups during the simulation to determine electron scattering after a collision.

[54] Each particle is then stepped through the “Null Collision Method” as summarized in Figure 6 to determine if it experienced a real collision in the time step, and if so, what type of collision occurred. If the electron experienced a real collision, the electron's energy and trajectory are updated corresponding to the type of collision as discussed in 9. After all electrons have been stepped through the collisional part of the model, each electron's parallel velocity (v_∥) is updated to reflect the acceleration due to the electric field E during the time step Δt using a first-order finite-difference representation of Newton's second law as v_∥^new = v_∥^old −EΔt.

[55] Diagnostic data is then saved to file to allow for time-dependent analysis of certain quantities to be calculated (16) after the simulation is complete (i.e., mean energy versus time, drift velocity versus time, rate coefficients, etc.). If a certain number of particles has been reached, the particle assembly is then remapped (14) in order to maintain reasonable computation times. It is then checked to determine if the simulation has converged to a steady state (i.e., mean energy and drift velocity remain constant over a certain time span). If the simulation has converged, then the simulation data is saved to file; if not, the simulation returns to the beginning of the next time step and the procedure is repeated until a converging solution is reached.

2.8. Model Diagnostics

[56] The electron energy distribution n() (normalized as n(ɛ) dɛ = 1) is obtained by sampling and averaging the electron assembly at several moments of time as the simulation reaches equilibrium. At each moment of time, the particles are divided into equally spaced bins along the energy ɛ axis and the particles weights W_j are summed to obtain the true number of electrons existing in each energy bin. The resulting functions are then normalized and averaged in time and the final distribution function n(ɛ) is normalized again

to ensure that n(ɛ) dɛ = 1, where ɛ_max is the maximum electron energy.

[57] The electron mean energy can be simply calculated by summing the energy of each particle included in the simulation

where W_j and ɛ_j are the weight and energy of the jth particle, respectively, N_t is the total number of particles in the simulation, and n_t is the total number of electrons which can be defined as

This calculation is performed at every time step and the result is used to determine if the simulation has reached an equilibrium state. Likewise, considering that the electron velocity is represented by its parallel and perpendicular velocity (Figure 9a, equations (38) through (40)), the drift velocity can be found by

The electron mobility can then be calculated as

[58] After a simulation has reached a steady state, the rate coefficients for various collision processes can be determined by simple counting procedures over a given time interval. First, consider a dummy variable C_i(t) which is used to count the occurrences of an ionization process over a given time interval. Each time an ionization collision occurs, the counter C_i(t) is incremented as

where W_j is the weight of the particle experiencing the collision. The rate coefficient for this ionization process can be obtained by selecting two moments in time t₁ and t₂ after the simulation has reached an equilibrium state. The ionization rate coefficient can then be calculated as

where n_t(t) is the total number of electrons as a function of time defined by equation (53). This calculation can be applied to any collision process (e.g., ν_a2 two-body attachment, B³Π_g excitation of N₂, C³Π_u excitation of N₂, etc.) for which a corresponding counter C_k(t) was maintained during the simulation. The meaning of the procedures defined by equations (56) and (55) can best be understood by direct integration of the dynamic equations describing the growth of the total number of electrons n_t due to ionization

or the growth of the total number of excited species n_k due to impact excitation by electrons

3.1. Comparisons With Previous Monte Carlo Model Results

[60] Figures 13, 14, and 15 show calculation comparisons of the electron energy distribution function n(ɛ) in N₂ gas at ground pressure from the current Monte Carlo model and the Monte Carlo model of Kunhardt and Tzeng [1986a]. For each simulation an assembly of 8000 electrons with an initially Maxwellian velocity distribution (electron temperature T_e = 0.5 eV) was placed under the influence of a uniform electric field and the simulation was performed until the assembly reached an equilibrium state. High-energy (ɛ > 500 eV) angular scattering of electrons was determined using equation (14) and remapping was performed when the total number of particles reached a value N_t > 15,000.

[61] The electron energy distribution function is shown for E/N values of 300 Td and 1500 Td (1 Td = 10⁻¹⁷ V cm²). Figures 13a and 14a are results presented by Kunhardt and Tzeng [1986a] corresponding to “model 3” of their study in which elastic and inelastic collisions were taken to be anisotropic and differential scattering cross sections were obtained from experiment and theory. For the same scattering consideration, the results of the present model are shown in Figures 13b and 14b. The energy distributions shown in Figures 13a and 13b demonstrate the same shape characteristics with the sharp peak in the number of electrons existing below 2 eV being ≈15% lower in the current calculations. This is most likely due to the use of updated collision cross section and differential cross section data in the current model, as compared to Kunhardt and Tzeng [1986a]. Similarly, Figures 14a and 14b also demonstrate the same overall shape with a ≈7% lower peak value for the same reason. It can be seen that for E/N = 1500 Td the high-energy tail of the distribution extends to energies ∼100 eV and with further inspection of the tail, Figures 15a and 15b show the population of thermal runaway electrons existing at energies ≫100 eV. The nonsmooth appearance of the distribution seen in Figure 15 is due to the small number of particles being sampled in these high-energy regions. Figure 15a corresponds to “Case 2” presented by Tzeng and Kunhardt [1986] in which the energy of secondary electrons emerging from ionizing collisions is determined from the differential ionization cross section of Opal et al. [1971] as discussed in 7. Kunhardt and Tzeng [1986a] emphasize that the treatment of angular scattering can have a great effect on the electron energy distribution function. While the variations displayed in the low-energy portion of the distribution are small, assumptions about angular scattering of electrons can significantly inhibit or promote the development of runaway as will be shown in 20, 22, and 23.

3.2. Comparisons With ELENDIF, Experimental Results, and Previous Models

[62] Figures 16, 17, and 18 show the energy distribution function n(ɛ), mean energy 〈ɛ〉, phase space velocity distribution, and drift velocity v_d of electrons in an air mixture at ground pressure consisting of 78.11% N₂, 20.91% O₂, and 0.98% Ar under the influence of applied electric fields E = 0.3E_k, E = 1.2E_k, and E = 20E_k, respectively, where E_k is the previously defined conventional breakdown field. As in 18, an assembly of 8000 electrons with an initially Maxwellian velocity distribution (electron temperature T_e = 0.5 eV) was used. Equation (14) was utilized to determine high-energy (ɛ > 500 eV for collisions with N₂ and ɛ > 1000 eV for collisions with O₂ and Ar) angular scattering of electrons and remapping was performed when the number of particles reached N_t > 15,000.

[63] The electron distribution function in Figures 16a, 17a, and 18a are compared to results obtained from the ELENDIF Boltzmann equation solver [Morgan and Penetrante, 1990]. It can be seen that for low field values E = 0.3E_k and E = 1.2E_k the energy distributions obtained from the current Monte Carlo model are in excellent agreement with results obtained from the ELENDIF software with differences of <5% and ≈5%, respectively. Also, Figures 16a and 17a provide an excellent insight into the dynamic friction force in air (2, Figure 2). Figure 16a shows that the majority of the electron population under the influence of an electric field E = 0.3E_k is maintained at energies <2 eV. Electrons are held at these low energies because of the N₂ vibrational processes with threshold energies ranging from 0.28 eV to 2.35 eV (see Table 1) constituting the “first bump” of the dynamic friction force (Figure 2). The electrons lose more energy to these vibrational collisions than they gain from the applied electric field. For an electric field E = 1.2E_k, however, it becomes possible for electrons to gain more energy from the electric field than they lose to these vibrational processes. As can be seen in Figure 17a, a large population of electrons remain at energies <2 eV due to vibrational collisions, but a significant number is able to penetrate through the vibrational barrier and is accelerated to energies >2 eV. These electrons, however, do not accelerate to higher energies because of the increased collision frequency and energy loss corresponding to the 100 eV “hump” of the dynamic friction force (Figure 2). Figure 18a shows the energy distribution of electrons when an electric field E = 20E_k, much greater than the thermal runaway field E_c, is applied and a large percentage of the total electron population is accelerated to energies >100 eV. It can also be seen in Figure 18a that the two-term spherical harmonic expansion [e.g., Raizer, 1991, chap. 5] used in the ELENDIF solution fails at this high electric field value and the Monte Carlo and ELENDIF results no longer agree. Figures 16b, 17b, and 18b show the simulation reaching an equilibrium state as the mean energy converges to a constant value for each of the electric field cases and Figures 16d, 17d, and 18d show the average drift velocity of electrons also converging. It can be seen from Figure 18b that the mean energy of electrons is only ∼60 eV, showing that even for an extremely high electric field E = 20E_k, it is still very difficult for electrons to accelerate to energies >100 eV and only a small fraction of them becomes runaway. Figures 16c, 17c, and 18c show the phase space velocity distributions of particles in the system. Figures 16c and 17c confirm that for electric fields E = 0.3E_k and E = 1.2E_k the velocity distributions remain largely isotropic, thus justifying the two-term spherical harmonic expansion used in the ELENDIF solution. On the contrary, Figure 18c shows that for E = 20E_k the velocity distribution becomes highly anisotropic because of thermal runaway electrons and therefore the two-term spherical harmonic expansion can no longer be utilized in the solution of the Boltzmann equation.

[64] Figures 19 through 22 show results of Monte Carlo model calculations compared to results obtained from ELENDIF, experimental swarm parameter measurements in dry air reported by Davies [1983], and from related analytical model approximations proposed by Aleksandrov et al. [1995] and Morrow and Lowke [1997] for a range of applied electric fields. Figure 19a shows that the electron mean energy calculated by the Monte Carlo model agrees well with the ELENDIF Boltzmann code at most electric fields, demonstrating only a slight deviation of ≈25% at E = 6 × 10⁷ V/m. The Monte Carlo calculation of electron mobility μ_e shown in Figure 19b shows excellent agreement with the experimental values of Davies [1983] and ELENDIF calculations for electric fields ranging from E = 5 × 10⁵ to E = 1 × 10⁷ V/m but demonstrates a percent difference of ≈38% with ELENDIF calculations at E = 6 × 10⁷ V/m. However, the Monte Carlo mobility calculation agrees well with the model approximation of Morrow and Lowke [1997] for electric fields >1×10⁷ V/m. It can be seen in Figure 20a that the ionization coefficient calculated by the Monte Carlo model demonstrates excellent agreement with experimental measurements [Davies, 1983], ELENDIF calculations, and the model of Aleksandrov et al. [1995] at all electric field values, but displays a difference of ≈68% with the model of Morrow and Lowke [1997] at E = 6 × 10⁷ V/m. The values calculated and reported for the two body attachment coefficient ν_a2 vary greatly, as can be seen in Figure 20b. The results of ν_a2 calculated by the Monte Carlo model agree well with the model of Aleksandrov et al. [1995] and ELENDIF calculations for electric fields >2 × 10⁶ V/m but are significantly higher at fields <2 × 10⁶ V/m. The Monte Carlo v_a2 calculations demonstrate a difference of up to ≈58% compared to experimental values, with the sharp decline in the ν_a2 values of Davies [1983] for E > 3 × 10⁶ V/m not being reproduced by Monte Carlo and ELENDIF solutions. This is most likely due to the presence of collisional detachment, which is not included in the Monte Carlo and ELENDIF calculations, impacting the experimental results [Davies, 1983]. It also can be seen from Figures 21 and 22 that the excitation coefficients of the B³Π_g and C³Π_u states of N₂ and the B²Σ_u⁺ state of N₂⁺ obtained from the Monte Carlo model agree well with results obtained from ELENDIF at high electric field values but show better agreement with the calculation of Aleksandrov et al. [1995] at lower electric fields. The B²Σ_u⁺ optical excitation coefficients shown in Figure 22 were calculated by obtaining the ionization rate coefficients of N₂ from the Monte Carlo model (see 16), ELENDIF calculations, and approximations of Aleksandrov et al. [1995], and then multiplying the resulting coefficients by a branching ratio (f_B = 0.145) presented by Van Zyl and Pendleton [1995].

[65] These results indicate that although the two-term spherical harmonic expansion is not generally valid at high electric field values, the above rate coefficients obtained from ELENDIF remain valid for electric fields up to 6 × 10⁷ V/m (≈20E_k). The differences between Monte Carlo and ELENDIF calculations of electron mean energy and mobility at high electric fields are most likely due to the electron distribution becoming anisotropic.

[66] We note that results of Figures 19 through 22 are presented in a reduced form using the ratio of the air density at ground pressure N₀ = 2.688 × 10¹⁹ cm⁻³ to the density at any other pressure N. This scaling is justified by the similarity properties of gas discharges [e.g., Liu and Pasko, 2004, and references therein] and can be directly used for calculation of physical quantities at any other air density N of interest.

[67] A MATLAB function, based on ELENDIF solutions, was used in the plotting of Figures 19 through 22. This function allows the calculation of many physical quantities and rate coefficients as a function of the reduced electric field in air and is freely available to readers as a supplement to this paper at http://pasko.ee.psu.edu/air.

3.3. Runaway Electron Development

[68] To examine the development of runaway electrons for various electric fields, simulations were conducted to perform direct comparisons with recent results presented by Bakhov et al. [2000]. An initial assembly of 1 eV electrons (N₀ = 5000) with velocities unidirectional with the electric force was placed in N₂ gas at ground pressure with applied electric fields of 400 kV/cm, 350 kV/cm, and 325 kV/cm, and simulated for times of t = 1 ns, t = 6 ns, and t = 25 ns, respectively. In accordance with Bakhov et al. [2000], the following nitrogen electronic states were considered (see Table 1): A³Σ_u⁺ (v = 5–9 and v = 10+), B³Π_g, W³Δ_u, B′³Σ_u⁻, C³Π_u, w¹Δ_u, a¹Π_g, and the sum of the remaining singlet states (threshold ɛ = 13 eV). Elastic and ionization collisions were also considered; however, while ionization was considered as an energy loss mechanism, secondary electrons were not added to the particle assembly. When electrons reached an energy equal to ɛ_max = 4000 eV, they were considered to be runaway and removed from the simulation.

[69] Figure 23 shows results of the Monte Carlo model of Bakhov et al. [2000], in which angular scattering of high-energy electrons after elastic and ionization collisions is determined from the electron data library (EEDL) [Perkins and Cullen, 1994] and is adopted to be forward after excitation collisions. Figure 24 shows similar results calculated by the Monte Carlo model discussed in this paper, where high-energy angular scattering has been determined by the differential scattering cross section approximation published by Kol'chuzhkin and Uchaikin [1978] (equation (14), Figure 24a) and the modified Rutherford cross section introduced in this paper (equation (16), Figure 24b). A direct inspection of the atomic nitrogen data presented in the EEDL library at selected energies of 1, 2.5, and 5 keV indicates that among the models discussed in 6 of this paper the corresponding differential cross sections are closest to the approximation given by equation (9).

[70] Recalling that low-energy (ɛ < 500 eV) electron scattering is determined by the experimentally determined differential cross sections listed in Table 5 for all simulations, the only difference between the simulations of Figure 24a and 24b is the differential scattering cross section used for energies >500 eV. It is obvious from the major differences shown in Figures 24a and 24b that the high-energy scattering of electrons is critical in the development of electron runaway. Comparison of Figures 4c and 4d show that the differential scattering cross section calculated from equation (16) possesses a much larger forward scattering characteristic than that of equation (14) for all electron energies. The effect of this forward scattering on the development of runaway electrons is well illustrated in Figure 24 and by comparing Figures 23 and 24 it can be seen that results obtained by Bakhov et al. [2000] using EEDL data are similar to those obtained with our model using the differential scattering cross section calculated from equation (16). It can also be seen from Figure 24a that assuming forward scattering after excitation collisions drastically increases the production of runaway electrons when equation (14) is used. This is not the case when the same assumption is used with equation (16) since equation (16) already possesses a strong forward scattering characteristic.

[71] Also of interest are the considerable differences in Figures 24a and 24b when forward scattering after excitation collisions is assumed. Since excitation collisions are forward for both cases, the elastic and ionization collisions must be responsible for the large differences in runaway electrons observed. While elastic collisions are the dominant collision process at low electron energies, the differences in Figures 24a and 24b are most likely due to the high-energy scattering of the primary electron after ionizing collisions. This leads to an interesting conclusion that the high-energy angular scattering from ionizing collisions is much more important in the development of electron runaway than high-energy scattering from excitation collisions.

4.1. Stationary Streamer Tip Model

[73] To demonstrate the importance of high-energy angular scattering in the development of runaway electrons in a streamer tip, four different cases utilizing the differential scattering cross section of Kol'chuzkin and Uchaikin [1978] (equation (14)) and the modified Rutherford cross section of equation (16) were considered. It can be seen from Figures 4 and 5 that these two cross sections possess very different scattering properties which greatly influence the development of runaway electrons as shown in Figure 24. As in 20, simulations were performed for both cross sections with and without the assumption of forward scattering after excitation collisions.

[74] The strong electric field associated with a streamer tip immediately preceding branching (Figure 1) is approximated by a 10 E_k square pulse spanning 1 m in z-space and is placed in various ambient electric fields E_amb (Figure 25a). We note that in accordance with similarity laws [e.g., Liu and Pasko, 2004], the streamer spatial scales (i.e., width of the streamer tip) and electric fields scale with air density N as ∼1/N and ∼N, respectively. Therefore the maximum energy, which can be gained by electrons in the streamer tip, remains the same for similar streamers at different air densities, or pressures (assuming constant gas temperature). The magnitude of the electric field in the streamer tip (10E_k = 2170 V/m) and the tip width (1 m) shown in Figure 25 correspond to a model streamer at 70 km altitude (Figure 1). As mentioned above, all results reported in this section can be easily generalized to streamers at other altitudes (pressures) of interest using similarity laws.

[75] The two quantities of interest are the runaway production rate ν_run [1/s] and the runaway flux Γ_run [1/m²/s]. The runaway production rate ν_run is defined as the rate at which electrons inside the streamer tip volume reach energies greater than a predefined runaway energy ɛ_run. The runaway flux Γ_run represents the number of electrons with energies greater than ɛ_run which pass through a unit area 1 m in front of the streamer tip per unit time. The distance of 1 m in front of the streamer tip was arbitrarily chosen such that the focus of our study is placed on the number of runaway electrons produced by streamer tips and not the future evolution of these electrons. Figure 26 shows the geometry and a sample calculation of ν_run and Γ_run for a case when the ambient field is assigned a value of E_amb = 4E_k and equation (16) was used to determine all angular scattering. Plotting the cumulative number of runaway electrons generated within the streamer tip volume N_run and the number of runaway electrons to pass through the unit area 1 m ahead of the streamer tip N_run-flux as functions of time, allows the rates ν′_run and Γ′_run to be calculated (see Figure 26). These rates may then be scaled to realistic values ν_run and Γ_run corresponding to a streamer tip electron density n_e at 70 km of 10¹² m⁻³ [Liu and Pasko, 2004] as

where n_tip is the electron density of the streamer tip in the simulation (see Figure 26). These values can then be generalized to streamers at other altitudes and pressures noting the similarity properties of streamers [Liu and Pasko, 2004, and references therein] following the relationship n_e = n_e0, where n_e0 = 2 × 10²⁰ m⁻³ is the reference streamer density at ground level, and N₀ and N are defined in the paragraph before last of 19.

[76] After performing a number of test simulations, it was determined that ν′_run and Γ′_run are independent of the 0.5 m E_amb field region behind the streamer tip (Figure 25a), thus allowing the simulation domain to be reduced to two electric field regions (Figure 25b). This simplification allows for more electrons to be placed in the 10 E_k streamer tip region, which improves the statistics of the ν′_run and Γ′_run calculations.

[77] In order to accurately represent the constant electron density of a propagating streamer tip [e.g., Liu and Pasko, 2004], the electron density was kept constant over the 10 E_k region throughout simulations using ten spatial bins (Figure 25c). After performing several test simulations, the electron density in the ambient field region n_amb was found to have no effect on ν_run and Γ_run and therefore was allowed to be variable. At the end of each time step, the density of electrons in each of the ten bins (n_bin) is calculated. The electron density in each bin is then compared with the desired constant n_tip to determine if there are too many or too few electrons in the bin. If n_bin > n_tip, then n_bin − n_tip electrons with energies less than 10 eV (ensuring that a runaway electron is not being repositioned) are randomly selected from the bin and 1 m is added to their current z-position thus moving them into the n_amb region. If n_bin < n_tip, then n_tip − n_bin electrons with energies less than 10 eV are randomly selected from the n_amb region and placed at the left boundary of the bin in question.

[78] This formulation results in a constant total number of electrons N_tip (see Figure 25c) being exposed to the high electric field of the streamer tip at all times throughout the simulation, but due to two-body and three-body attachment collisions with O₂ (see Table 2) the total number of electrons in the low field region N_amb decreases slightly over time. Initial values of n_tip are normally chosen to be 150,000, 100,000, 75,000, 50,000, or 25,000 m⁻³ depending on the resolution needed to obtain satisfactory statistics for various electric field and scattering cases, and n_amb is initially 20,000 m⁻³ for all simulations.

[79] It should be noted that while ionization is considered as an energy loss mechanism for primary electrons, secondary electrons are not added to the simulation of the stationary streamer tip, and the one-dimensional remapping of particles discussed in 14 was not utilized in the simulations. For an electric field of 10E_k the ionization coefficient is ∼7 × 10¹¹ s⁻¹ (see Figure 20), which leads to an enormous multiplication of electrons. If these secondary electrons are added to the simulation, as with the zero-dimensional model, the particle weights and electron density grow at an extremely fast pace. This fast growth is believed to be an inaccurate representation of the actual physical system in which the electron density in the streamer tip is dynamically maintained at relatively stable levels. The focus of this study is therefore placed on determining the number of initially cold electrons out of a known constant electron density to reach runaway energies. This evolution of individually tracked runaway electrons can therefore be treated without consideration of the avalanche multiplication of electrons due to ionization. This assumption, however, relies heavily on the fact that most secondary electrons emerge from ionizing collisions with low energies. Simulations have shown that a small percentage of secondaries do possess energies >200 eV, thus placing them over the “hump” of the dynamic friction force (Figure 2). These high-energy secondaries entering an electric field ∼10 E_k have an increased probability of becoming runaway electrons and would impact the values of ν_run and Γ_run. Several simulations were performed in which secondary electrons emerging from ionizing collisions with energies >200 eV were added to the simulation and values of Γ_run were observed to increase by ≈10%.

[80] The energy required for an electron to be considered runaway ɛ_run is directly determined from the dynamic friction force of electrons in air (Figure 2). Figure 27a shows the values of ɛ_run corresponding to ambient electric fields of 1 E_k, 2 E_k, 3 E_k, 4 E_k, and 5 E_k. It can be seen from Figure 2a that these ambient fields alone are not strong enough to overcome the friction force F_D and accelerate low-energy electrons to energies >100 eV. An electric field of 10E_k, however, falls above the peak of F_D and it becomes possible for electrons exposed to a 10 E_k field to be accelerated to energies >100 eV. In view of the electric field configuration shown in Figure 25b, the 10E_k streamer tip field extending over 1 m length is solely responsible for all runaway electrons produced during a simulation, and the magnitude of the ambient field E_amb defines the value ɛ_run. To better understand this relationship, consider a case with E_amb = 3 E_k depicted in Figure 27b. An electron emerging from the 10 E_k streamer tip region into the ambient field region may be considered to fall into one of three categories: the electron has an energy ɛ < ɛ_run, ɛ ∼ ɛ_run, or ɛ > ɛ_run. If the electron has not gained sufficient energy from the 10 E_k streamer tip (ɛ < ɛ_run), the electron's motion will be dominated by collisions and the friction force will cause the electron to decelerate to low energies. If the electron has gained enough energy from the streamer tip such that ɛ > ɛ_run, the probability of the electron colliding with neutral particles has significantly decreased (F_D < 3 E_k) and the electron will continue to be accelerated to higher energies by the 3 E_k ambient field. If the electron has gained an energy ɛ ∼ ɛ_run, the electron may be accelerated to higher energies by the ambient field or return to low energies depending on the collisions and angular scattering it experiences.

[81] Table 7 lists the values of ɛ_run for various ambient electric fields and the corresponding electron velocities. It should be noted that the maximum energy which an electron experiencing no collisions could gain from the 1 m long streamer tip (see Figure 26) with 10 E_k field at 70 km is 2170 eV. We emphasize again that this energy will be exactly the same for similar streamers at any pressure/altitude of interest. This energy corresponds to an ambient field of 2.32 E_k. Therefore it is impossible to achieve runaway electrons for ambient fields <2.32 E_k in the simulations. In reality, when collisions are considered, electrons only gain a fraction of the total 2170 eV and no runaway electrons were observed in any simulations with ambient electric fields 2.5 E_k. Simulations were therefore performed for ambient field values ranging from 2.5 E_k to 5 E_k as shown schematically in Figure 25b.

[82] We emphasize that changes in the width of the model streamer and the magnitude of the tip electric field would affect the numerical values of ambient fields at which runaway electrons can be sustained, and these values can be easily evaluated using procedures outlined above.

[83] Figures 28 and 29 show simulation results for the runaway production rate ν_run and the runaway flux out of the streamer tip Γ_run for various E_amb values (E_tip = 10 E_k for all simulations). Figure 28a shows results obtained when high-energy angular scattering after all collisions was determined using equation (16), and Figure 28b shows the same case except high-energy scattering after excitation collisions was considered to be forward (χ = 0). Figures 29a and 29b show results for the same scenarios, respectively, utilizing equation (14) to determine high-energy scattering instead of equation (16). The solid lines in both Figures 28 and 29 were obtained using a spline interpolation between Monte Carlo data points.

[84] The variations in ν_run and Γ_run at different ambient fields can be largely attributed to the runaway energy ɛ_run being variable with the ambient field (Table 7). For example, consider ambient fields of 5E_k and 3 E_k corresponding to ɛ_run = 630 eV and ɛ_run = 1515 eV, respectively. A much larger fraction of the electron distribution exists at energies >630 eV than at energies >1515 eV, thus explaining the notable increase in ν_run and Γ_run with the increase in ambient field. We emphasize again that the generation of runaway electrons is entirely due to the 10 E_k streamer tip field. If ɛ_run were assigned a constant value for all ambient field cases, ν_run would also remain constant for all ambient fields since the 10 E_k streamer tip remains unchanged.

[85] It should be noted that not all electrons which achieve energies >ɛ_run can be considered true runaway electrons. The direction in which the electrons are traveling with respect to the electric field also plays a major role in runaway development. Electrons traveling in a direction opposing the applied electric field gain additional energy from the electric field, whereas electrons traveling in the same direction as the electric field lose energy. Therefore it is possible for electrons with energies >ɛ_run to have v_∥ velocity components in the same direction as the electric field and subsequently be losing energy. Clearly, these electrons do not fit the definition of runaway electrons since they are in fact losing energy to the electric field. This illustrates the reason why the angular scattering of electrons becomes so important in the development of runaway electrons. When electrons are subject to significant angular scattering, they are more prone to obtain a v_∥ which is in the same direction as the applied electric field. Adversely, if electrons scatter primarily forward after collisions (in the direction of electron drift and opposing the direction of the electric field), it becomes easier for them to become runaway electrons because they will always be gaining energy from the applied electric field. In order to more accurately differentiate between true runaway electrons which are traveling in the direction opposite to the electric field and high-energy electrons which are traveling in the same direction as the electric field, a filtering procedure must be introduced into the simulations. Figure 30 shows the calculation of two fluxes of high-energy electrons obtained from a model simulation, Γ₀ which is leaving the streamer tip at z = 1 m and Γ_run which was defined previously. The results shown in Figure 30 are for a case when E_amb = 4E_k, all angular scattering was determined using equation (16), and N_tip = 50,000. The difference between Γ₀ and Γ_run arises due to the 1 m ambient field region ahead of the streamer tip and the deceleration of a fraction of high-energy electrons within it. As discussed above, an electron exiting the streamer tip with an energy ɛ ∼ ɛ_run may either continue to be accelerated to high energies by the ambient electric field or be returned to low energies due to collisions possessing strong backscatter components. Therefore the ambient field region acts as a filter separating the two groups of electrons. An electron which still possesses an energy greater than ɛ_run after traveling 1 m in the ambient field region continues to gain more energy from the ambient field than it lost to collisions, consequently obeying the definition of a runaway electron as presented in 2. It can be seen from Figure 30 that for this case only 39% of the electrons which achieved energies >ɛ_run at the streamer tip (z = 1 m), still possessed energies >ɛ_run after traveling through the 1 m ambient field region, the other 61% of electrons being decelerated by collisions in the ambient field region.

[86] Generally, the difference between Γ₀ and Γ_run grows larger as the ambient field intensity increases (ɛ_run decreases). At low ambient fields (high ɛ_run) electrons which achieve energies >ɛ_run and contribute to Γ₀ are at high enough energies such that they scatter predominantly forward after collisions. Therefore a large fraction of electrons achieving these high energies will continue to scatter forward in the ambient field region and accelerate to higher energies, thus contributing to Γ_run. Adversely, for high ambient fields (low ɛ_run) electrons are counted as runaway within the streamer tip at much lower energies and do not possess the same forward scattering properties. Therefore a large number of these electrons will be subject to significant angular scattering upon entering the ambient field region and will decelerate to low energies bringing about a larger difference in Γ₀ and Γ_run.

[87] Table 8 lists the average energies of electrons which were counted in the calculation of Γ_run and the maximum electron energy observed of those electrons. It becomes clear by comparing the values listed in Table 8 with the runaway energies ɛ_run listed in Table 7 that a portion of the electrons which were accelerated to energies >ɛ_run by the streamer tip continued to gain energy from the ambient electric field outside of the streamer tip, thus satisfying the definition of a true runaway electron. For most electric field cases, electron energies >2 keV were observed at the left boundary 1 m ahead of the streamer tip. However, it should be noted that while many electrons counted in the calculation of Γ_run sustained energies significantly greater than ɛ_run, some electrons still possessed energies ∼ɛ_run and therefore may be decelerated to low energies by collisions at distances greater than 1 m ahead of the streamer tip. The average values listed in Table 8 are dependent on the number of electrons included in the simulation. Since various electron densities were used for different ambient field and scattering cases (i.e., statistics are generally better for higher ambient fields), the accuracy of the values listed in the tables may vary with the ambient field.

[88] One of the most interesting aspects of Figures 28 and 29 is the existence of a cutoff electric field at which ν_run and Γ_run rapidly drop to zero. In Figure 28a, in which equation (16) was used to determine all high-energy angular scattering, the cutoff fields are E_amb ≈ 2.5E_k for ν_run and E_amb ≈ 2.55E_k for Γ_run. When excitation collisions are assumed to be forward for the same case, the cutoff electric field for both ν_run and Γ_run is E_amb ≈ 2.5E_k as can be seen in Figure 28b. This would suggest that the assumption of forward scattering after excitation collisions does not affect the cutoff electric field of the runaway production rate ν_run but does impact the cutoff field of the runaway flux Γ_run. It can also be seen from Figures 28a and 29a that at high electric fields (>E_amb = 3E_k), the forward scattering assumption has a minimal effect on the magnitude of ν_run and Γ_run, possessing a maximum percent difference of ≈14% and ≈19%, respectively. However, for electric fields E_amb = 3E_k, the percent differences between Figures 28a and 29b are ≈65% for ν_run and ≈137% for Γ_run, demonstrating that angular scattering after excitation collisions plays a larger role when lower ambient fields (higher ɛ_run) are applied.

[89] Similarly, Figure 29a, in which equation (14) was used to determine all high-energy angular scattering, exhibits cutoff electric fields of E_amb ≈ 2.6E_k and E_amb ≈3E_k for ν_run and Γ_run, respectively. When scattering after excitation collisions is assumed to be forward, the cutoff field of Γ_run decreases from E_amb ≈ 3E_k to E_amb ≈2.6E_k while the cutoff field of ν_run remains E_amb ≈ 2.6E_k. Similar to Figure 28, the differences between the ν_run values of Figures 29a and 29b are relatively small, with a maximum difference of ≈83%. On the contrary, the difference between the Γ_run values plotted in Figures 29a and 29b is dramatic, possessing percent differences >100% at all electric field values. The fact that differences between Figures 29a and 29b are much greater than the differences between Figures 28a and 28b demonstrates that assuming forward scattering after excitation collisions has a much greater effect on runaway electrons when equation (14) is used to calculate high-energy angular scattering and corresponds well with the results reported in Figure 24.

[90] As discussed in 6 and shown in Figure 24, the differential scattering cross section used to describe high-energy electron scattering can have a major impact on the development of runaway electrons. This can further be seen by comparing the results of Figures 28a to 29a and the results of Figures 28b to 29b. As discussed above, the most notable difference is the variation in the cutoff electric field at which ν_run and Γ_run rapidly drop to zero. In addition to the variation in the cutoff electric field, the magnitudes of ν_run and Γ_run are noticeably greater when the differential cross section of equation (16) is used. Comparing Figures 28a to 29a, at ambient electric fields ≥3.5E_k, the values of Γ_run in Figure 28a are ≈25 to 45 times greater than those of Figure 29a. When scattering after excitation collisions is assumed to be forward, the Γ_run values of Figure 28b at ambient fields ≥2.8E_k are ≈8 to 16 times those of Figure 29b. The differences in ν_run are not as drastic, with a maximum percent difference between the values shown in Figure 28a and Figure 29a of ≈150% at E_amb = 2.8 E_k decreasing to ≈40% at E_amb = 5 E_k. Similarly, comparing Figures 28b to 29b results in percent differences of ≈140% at E_amb = 2.8 E_k decreasing to ≈25% at E_amb = 5 E_k. Since the differences in ν_run are not as great as those in Γ_run, one may conclude that the angular scattering of electrons once they enter the ambient field region ahead of the streamer tip is much more critical for sustaining runaway electrons than the scattering of electrons within the high field streamer tip.

[91] The difference in ν_run between the two cross section cases (equations (14) and (16)) decreases at higher electric field due to the fact that as E_amb increases, ɛ_run decreases. For example, when E_amb = 5 E_k, ɛ_run = 630 eV, so recalling from 5 and 6 that electron scattering from collisions with N₂ at energies ɛ < 500 eV and from collisions with O₂ and Ar at energies ɛ < 1000 eV are treated using experimental cross section data for all simulations, the scattering of electrons which achieve energies >ɛ_run has only been determined using equations (14) and (16) over the small energy range of 500 to 630 eV. On the contrary, when E_amb = 3 E_k (ɛ_run = 1515eV), scattering of electrons reaching energies >ɛ_run have been treated using equations (14) and (16) over much larger energy ranges of 500 to 1515 eV for collisions with N₂ and 1000 to 1515 eV for collisions with O₂ and Ar, therefore the effect of high-energy scattering approximations is much greater for lower electric fields (higher ɛ_run).

[92] Figure 31 shows the electron energy distribution function, drift velocity, and mean energy as well as an electron position versus electron energy phase space plot for six different regions in the one-dimensional simulation. The ambient electric field was chosen to be E_amb = 4E_k and equation (16) was used to calculate high-energy electron scattering. It can be seen in Figure 31 that the electron population within the streamer tip maintains a drift velocity ∼10⁶ m/s and a mean energy ∼30 eV. Just ahead of the streamer tip, the drift velocity and mean energy were found to be ≈5.24 × 10⁵ m/s and ≈17.0 eV, respectively, demonstrating a 62% and 57% decrease from just inside the streamer tip. These values continue to decrease to values of ≈3.94 × 10⁵ m/s and ≈10.6 eV at a distance ≈1 m ahead of the streamer tip. From the electron position versus energy phase space plots, it can be seen that while the bulk of the electron population remains at low energies, a small number of electrons are accelerated to runaway energies by the streamer tip and continue to runaway as they travel through the ambient field region. It can also be seen from the phase space plots that a number of electrons exiting the streamer tip with high energies eventually decelerate as they move further away from the streamer tip.

4.2. Self-Acceleration of Electrons

[93] Simulation results in 22 have been presented for a case when the streamer tip was assumed to be stationary for computational simplicity. In reality, streamers can propagate at velocities on the order of one tenth of the speed of light. Stanley et al. [1999], Moudry et al. [2002, 2003], and McHarg et al. [2002] have reported streamer speeds in sprites of ∼1.2 × 10⁷, ∼3 × 10⁷, and ∼5.3 × 10⁷ m/s, respectively. The electron energies corresponding to these speeds are 409 eV, 2559 eV, and 7985 eV (see Table 7). Therefore for an electron to exit the streamer tip ahead of the propagating streamer it must achieve a velocity greater than that of the streamer's propagation. In view of the simulation setup presented in 22, since the maximum energy an electron can gain in the 1 m streamer tip region is 2170 eV, it would seem to be impossible for an electron to accelerate ahead of streamers propagating at the latter two speeds. However, a situation may arise during the propagation of a negative streamer in which a fraction of the total electron population travels with the propagating streamer tip. This is known as self-acceleration of electrons and was first proposed by Babich [1982]. If this were to occur, electrons would be exposed to the high streamer tip field for an extended amount of time and may gain enough additional energy from the high electric field to accelerate ahead of the streamer tip and may make it possible for runaway electrons to be observed for ambient electric fields <2.5E_k.

[94] In order to study the effects of self-acceleration of electrons, in this section the simulation setup presented in 22 is modified to account for the motion of the streamer tip. Using the same electric field configuration as discussed in 22 (see Figure 25b), the relative velocity of the streamer tip may be added to the simulation by introducing an extra term to equation (60) as

where v_str is the velocity of the streamer. Therefore in each time step the electrons will travel a distance v_∥ · Δt according to their parallel velocity component but will also be moved a distance v_str · Δt back in the simulation space to reflect the movement of the streamer. Essentially, the electrons are moving at their own speed, and the simulation domain is moving at the speed of the streamer v_str.

[95] The electron density is initially constant over the entire 2 m span and remains so throughout the simulation without any remapping or repositioning of the electrons being necessary. When runaway electrons are lost from the simulation through the 1 m² area ahead the streamer tip, low-energy electrons are randomly positioned back into the simulation to replace them. When low-energy electrons pass through the z = 0 m left side of the simulation domain (see Figure 25b), they are repositioned at z = 2 m and their energies are reset to the mean energy of electrons corresponding to the E_amb electric field. Similarly, electrons passing through the z = 2 m boundary are repositioned at z = 0 m and their energies are reset to the mean energy of electrons corresponding to the E_tip electric field. The simulation may simply be viewed as a 1 m long electric field pulse traveling at a speed v_str through a constant density of electrons which are traveling at a drift velocity v_d determined by the applied ambient field E_amb. Since v_d is generally much less than v_str, the bulk of the electron population will simply be passed over by the streamer tip and repositioned at z = 2 m as described above. However, depending on the streamer speed v_str, it is possible for a small fraction of the total electron population to actually achieve an energy necessary to travel with the streamer tip and experience the self-acceleration phenomenon. As with the stationary streamer tip model, secondary electrons emerging from ionizing collisions were not added to the simulation. If secondary electrons were included, similar increases of Γ_run could be expected as discussed in the previous section.

[96] Simulations were first performed for streamer tip fields of 10E_k and streamer speeds of v_str = 3 × 10⁷ and v_str = 5.3 × 10⁷ m/s. For these streamer speeds no runaway electrons were observed as the streamer propagated too fast to allow adequate time for electrons to gain sufficient energy from the 10E_k electric field and travel with the streamer tip. The electric field pulse essentially flew right by the background electron density. However, when the streamer speed was set to v_str = 1.2 × 10⁷ m/s, the results were drastically different.

[97] Figure 32 shows the runaway flux Γ_run for various ambient electric fields (for all simulations E_tip = 10 E_k and equation (16) was used to calculate angular scattering of high-energy electrons after all collisions). As can be seen from Figure 32, the ambient fields at which 10 E_k streamer tips can produce thermal runaway electrons extend all the way down to the conventional breakdown field E_k. Considering the major differences between the values of Γ_run shown in Figure 28a and Figure 32, it can easily be concluded that the self-acceleration of electrons can play a major role in the production of thermal runaway electrons by streamer tips. In addition to the increased fluxes of thermal runaway electrons, the energies achieved by the runaway electrons also increased significantly as can be seen in Table 9. Once again it should be noted that the variation in flux for various ambient fields arises due to the runaway energy ɛ_run being defined by the ambient field as discussed in 22. The average runaway electron energies increase for electric fields <2E_k due to ɛ_run increasing (see Table 7), therefore eliminating lower energy runaway electrons from the averaging which would have been included for higher ambient fields (lower ɛ_run). The increase in ɛ_max for electric fields <2E_k can be attributed to increased electron densities used in the simulations which in turn increase the chance that one electron out of the assembly may be accelerated to very high energies.

[98] In addition to the simulations shown in Figure 32, simulations were also performed for a streamer tip field of E_tip = 8E_k and streamer speed of v_str = 1.2 × 10⁷ m/s. Using the same procedures as outlined above and equation (16) to describe high-energy electron scattering, runaway fluxes of Γ_run = 6.12 × 10¹² 1/m²/s and Γ_run = 1.17 × 10¹³ 1/m²/s were obtained for ambient fields of E_amb = 2E_k and E_amb = 3E_k, respectively. The average and maximum energies were found to be ɛ_avg = 3707 eV and ɛ_max = 4681 eV for E_amb = 2E_k and ɛ_avg = 3909 eV and ɛ_max = 4228 eV for E_amb = 2E_k. For this E_tip = 8E_k case, simulations with ambient electric fields <2E_k produced no runaway electrons.

[99] The results presented above suggest that there exists a range of streamer propagation speeds and electric field configurations capable of producing a significant number of thermal runaway electrons. Performing simple analytical calculations (with electron-neutral collisions neglected), it can be found that for a 1 m in width streamer tip of 10E_k there exist a cutoff streamer speed v_str at which electrons cannot experience self-accelerating effects. It can be shown for this streamer configuration that at streamer speeds v_str > 2.75 × 10⁷ m/s no runaway electrons will flux through the tip of streamer since this speed does not allow adequate time for low-energy electrons to gain an energy substantial enough to travel with the propagating streamer. Conversely, self-accelerating effects can be observed for a wide range of lower streamer propagation speeds; however, the maximum energy an electron can gain under the most favorable conditions never exceeds ∼9 keV. For v_str = 1.2 × 10⁷ m/s (corresponding to the results discussed above), it can be found that the maximum possible energy of an electron leaving the streamer tip after experiencing self-acceleration is ∼4.5 keV. This value corresponds well with most of the average runaway energy ɛ_avg values listed in Table 9. It should be noted that the value of 4.5 keV estimated above from an analytical formulation is calculated for electrons exiting the streamer tip at z = 1 m and does not account for electrons' perpendicular velocity components v_⊥. Electrons which were observed in simulations with energies ∼8 keV obtained perpendicular velocities v_⊥ > v_∥. These electrons most likely moved in and out of the streamer tip region several times before finally exiting through the z = 2 m boundary, allowing them to gain additional energy.

4.3. Fluxes of Thermal Runaway Electrons in Lightning Leaders

[100] The results for ν_run and Γ_run presented in 22 and 23 are for altitudes of 70 km corresponding to sprites. However, these values may also be scaled using the similarity properties of streamers [e.g., Liu and Pasko, 2004, and references therein] to other altitudes of interest such as ground level [Dwyer et al., 2005] and 16 km [Dwyer, 2005b], which may correspond to the leader streamer zone of lightning. Assuming the neutral atmospheric densities N and streamer densities n_e listed in Table 10 for altitudes of 0, 16, and 70 km, the values of ν_run and Γ_run reported in Figures 28, 29, and 30 may be scaled to 0 and 16 km simply by multiplying the values by scaling factors and , respectively.

[101] As discussed in 1, a leader fuses the currents of numerous streamers that start from the tip, into a single channel [Raizer, 1991, p. 367]. Now, consider a leader with a typical current ∼100 A [Raizer, 1991, p. 372] propagating at an altitude of 16 km. The current of a typical streamer may be approximated by

where v_d is the drift velocity of electrons (v_d = μ_e∣E∣) and r_s is the radius of the streamer. For a streamer radius r_s = 2.3 × 10⁻³ m and the streamer channel field ∼ E_k, the typical streamer current would be I_s ≈ 1.8 A. Dividing the leader current of 100 A by the individual streamer current of 1.8 A would result in an estimation of 55 streamers ahead of the leader tip. We note that 55 is a very conservative estimate. It is known that up to 10⁵ streamers can exist at any moment of time in the leader streamer zone during the quasi-stationary stages of leader propagation [Bazelyan and Raizer, 2000, p. 70]. The assumed number (55) should adequately represent the total number of streamers attached to the leader head during a transient (∼1 μs) negative corona flash stage of the negative leader development discussed in 26.

[102] Scaling the value of Γ_run from Figure 32 at E_amb = 1.5E_k to 16 km results in Γ_run = 1.44 × 10²¹ 1/m²/s and multiplying by 55 streamers gives a total runaway flux of ≈7.92 × 10²² 1/m²/s associated with the leader streamer zone at 16 km. Considering that the cross-sectional area of a streamer with radius r_s = 2.3×10⁻³ m is 1.66 × 10⁻⁵ m² and multiplying the cross-sectional area of one streamer by 55 total streamers results in a total area of 9.14 × 10⁻⁴ m². The runaway flux Γ_run = 1.44 × 10²¹ 1/m²/s may then be multiplied by the total cross-sectional area 9.14 × 10⁻⁴ m² of the 55 streamers resulting in a total of ≈10¹⁸ runaway electrons emitted by the leader head per second. Noting that the streamer electron number density and the flux Γ_run scale with the atmospheric neutral density as ∼N² and the streamer cross-sectional area scales as ∼1/N² (since r_s ∼ 1/N), the estimate of the total number of runaway electrons produced by the leader tip per second presented above (i.e., 10¹⁸ s⁻¹) is valid for a 100 A leader at any other air pressure/altitude of interest.

[103] We emphasize that the calculations above are based on runaway fluxes calculated at a fixed distance ahead of the streamer tips, and at other distances the number of runaway electrons may vary greatly from the values reported in this paper depending on particular details of geometry and time dynamics of the driving electric field. These estimates nevertheless demonstrate that a significant number of thermal runaway electrons may be generated by the leader streamer zone of lightning, which could contribute to the generation of recently reported X-ray and gamma ray bursts observed in association with lightning discharges [Moore et al., 2001; Dwyer et al., 2003, 2004a, 2004b, 2005; Smith et al., 2005].

[104] Concerning the theory of runaway breakdown, the future dynamics of these electrons at greater distances from the streamer tip must be considered. Depending on the geometry and time dynamics of the electric field ahead of the streamer tip, runaway electrons with energies ∼2 keV generated by the streamer tip may either decelerate to low energies or continue to gain energy ∼MeV. These MeV electrons may then contribute to relativistic runaway avalanches, as described by Gurevich et al. [1992], Roussel-Dupre et al. [1994], Lehtinen et al. [1999], Gurevich and Zybin [2001], Dwyer [2003], and references cited therein (see discussion in 26). A more complex model should be introduced to link thermal runaway electrons generated by streamer tips to the relativistic runaway electron avalanche model proposed by Gurevich et al. [1992] and determine the fraction of thermal runaway electrons which obtain energies ∼MeV for various electric field configurations.